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Abstract 

' A strong variable gravitational field of the very early Universe inevitably generates relic gravitational waves 

_ by amplifying their zero-point quantum oscillations. We begin our discussion by contrasting the concepts 

O I of relic gravitational waves and inflationary 'tensor modes'. We explain and summarize the properties 

^ ■ of relic gravitational waves that are needed to derive their effects on CMB temperature and polarization 

anisotropics. The radiation field is characterized by four invariants /, V, E, B. We reduce the radiative 
transfer equations to a single integral equation of Voltairre type and solve it analytically and numerically. 
, We formulate the correlation functions Cf^ for X, X' ~ T,E,B and derive their amplitudes, shapes and 



oscillatory features. Although all of our main conclusions are supported by exact numerical calculations, we 
obtain them, in effect, analytically by developing and using accurate approximations. We show that the TE 
correlation at lower i"s must be negative (i.e. an anticorrelation) , if it is caused by gravitational waves, and 
positive if it is caused by density perturbations. This difference in TE correlation may be a signature more 
valuable observationally than the lack or presence of the BE correlation, since the TE signal is about 100 
times stronger than the expected BB signal. We discuss the detection by WMAP of the TE anticorrelation 
at £ ~ 30 and show that such an anticorrelation is possible only in the presence of a significant amount 
of relic gravitational waves (within the framework of all other common assumptions). We propose models 
containing considerable amounts of relic gravitational waves that are consistent with the measured TT, TE 
and EE correlations. 

PACS numbers: 98.70.Vc, 98.80.Cq, 04.30.-w 
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I. INTRODUCTION 



The detection of primordial gravitational waves is rightly considered a highest priority task for 
the upcoming observational missions 1]. Relic gravitational waves are inevitably generated by 
strong variable gravitational field of the very early Universe. The generating mechanism is the 

n 

superadiabatic (parametric) amplification of the waves' zero-point quantum oscillations ;2]. In 
contrast to other known massless particles, the coupling of gravitational waves to the 'external' 
gravitational field is such that they could be amplified or generated in a homogeneous isotropic 
FLRW universe. This conclusion, at the time of its formulation, was on a collision course with the 
dominating theoretical doctrine. At that time, it was believed that the gravitational waves could 
not be generated in a FLRW universe, and the possibility of their generation required the early 
Universe to be strongly anisotropic (see, for example, P]). 

The generating mechanism itself relies only on the validity of general relativity and quantum 
mechanics. But the amount and spectral content of relic gravitational waves depend on a specific 
evolution of the cosmological scale factor (classical 'pumping' gravitational field) a{rj). The theory 
was applied to a variety of a(rj), including those that are now called infiationary models (for a 
sample of possible spectra of relic gravitational waves, see Fig. 4 in Ref. 4]). If a unique a{r]) 
were known in advance from some fundamental 'theory of everything', we would have derived the 
properties of the today's signal with no ambiguity. In the absence of such a theory, we have to 
use the available partial information in order to reduce the number of options. This allows us to 
evaluate the level of the expected signals in various frequency intervals. The prize is very high - 
the actual detection of a particular background of relic gravitational waves will provide us with 
the unique clue to the 'birth' of the Universe and its very early dynamical behaviour. 

A crucial assumption that we make in this and previous studies is that the observed large- 
scale cosmic microwave background (CMB) anisotropies are caused by cosmological perturbations 
of quantum-mechanical origin. If this is true, then general relativity and quantum mechanics 
tell us that relic gravitational waves should be a significant, if not a dominant, contributor to 
the observed large-scale anisotropies. From the existing data on the amplitude and spectrum of 
the CMB fiuctuations we infer the amplitude and spectral slope of the very long relic gravitational 
waves. We then derive detailed predictions for indirect and direct observations of relic gravitational 
waves in various frequency bands. 

At this point, it is important to clarify the difference between the concepts of relic gravitational 
waves and what is now called inflationary gravitational waves. The statements about infiationary 
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gravitational waves ('tensor modes') are based on the inflation theory. This theory assumes that 
the evolution of the very early Universe was driven by a scalar field, coupled to gravity in a special 
manner. The theory does not deny the correctness of the previously performed calculations for 
relic gravitational waves. However, the infiationary theory proposes its own way of calculating the 
generation of density perturbations ('scalar modes'). The infiationary theory appeals exactly to 
the same mechanism of superadiabatic (parametric) amplification of quantum vacuum fluctuations, 
that is responsible for the generation of relic gravitational waves, but enforces very peculiar initial 
conditions in the 'scalar modes' calculations. 

According to the inflationary initial conditions, the amplitudes of the 'gauge-invariant' metric 
perturbations (" associated with the 'scalar modes' (or, in other words, the amplitudes of the 
curvature perturbations called ( or TZ) can be arbitrarily large from the very beginning of their 
evolution. Moreover, the theory demands that these amplitudes must be infinitely large in the 
limit of the deSitter expansion law a(r/) oc which is responsible for the generation of a fiat 

(Harrison-Zeldovich-Peebles, 'scale-invariant') primordial spectrum, with the spectral index n = 1. 
At the same time, the amplitudes of the generated gravitational waves are finite and small for 
all spectral indices, including n = 1. Since both, gravitational waves and density perturbations, 
produce CMB anisotropies and we see them small today, the infiationary theory substitutes (for 
'consistency') its prediction of infinitely large amplitudes of density perturbations, in the limit 
n = 1, by the claim that it is the amount of primordial gravitational waves, expressed in terms 
of the 'tensor /scalar ratio r', that should be zero, r = 0. For a detailed critical analysis of the 
inflationary conclusions, see P]; for arguments aimed at defending those conclusions, see Q]. 

The science motivations and CMB data analysis pipelines, designed to evaluate the gravitational 
wave contribution, are usually based on inflationary formulas ^jQjyj^^l- In particular, according 
to the inflation theory, the primordial power spectrum of density perturbations has the form (it 
follows, for example, from Eqs. (2.12a) and (2.12b) in Ref.|lo| or from Eqs. (18) and (19) (or 
(A12) and (A13)) in Ref.Q): 

Despite the fact that this spectrum diverges at r = (n^ = 1, nj = 0) the CMB data analysts 
persistently claim that the inflation theory is in spectacular agreement with observations and the 
CMB data are perfectly well consistent with r = (the published confldence level contours always 
include r = and are typically centered at that point). 

Our analysis in this paper, based on general relativity and quantum mechanics, is aimed at 
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showing that there is evidence of signatures of rehc gravitational waves in the ah'eady available 
CMB data. We also make predictions for some future experiments and observations. 

The plan of the paper is as follows. In Sec. ^ we summarize the properties of a random 
background of relic gravitational waves that are needed for CMB calculations. The emphasis is on 
the gravitational wave (g.w.) mode functions, power spectra, and statistical relations. In Sec. II 111 
we discuss the general equations of radiative transfer and explain the existence of four invariants 
/, V, E, B that fully characterize the radiation field. We formulate the linearized equations in the 
presence of a single Fourier mode of gravitational waves. We prove that there exists a choice 
of variables that reduces the problem of temperature and polarization anisotropics to only two 
functions of time a(r/) and (3{ri). 

Section HVl is devoted to further analysis of the radiative transfer equations. The main result of 
this section is the reduction of coupled integro-differential equations to a single integral equation of 
Voltairre type. Essentially, the entire problem of the CMB polarization is reduced to a single inte- 
gral equation. This allows us to use simple analytical approximations and give transparent physical 
interpretation. In Sec. we generalize the analysis to a superposition of random Fourier modes 
with arbitrary wavevectors. We derive (and partially rederive the previously known) expressions for 
multipole coefficients af^ of the radiation field. We show that the statistical properties of the mul- 
tipole coefficients are fully determined by the statistical properties of the underlying gravitational 
perturbations. This section contains the expressions for general correlation and cross-correlation 
functions Cf-^' for invariants of the radiation field. 

We work out astrophysical applications in Sec. IVII and Sec. IVIII where we discuss the effects of 
recombination and reionization era, respectively. Although all our main conclusions are supported 
by exact numerical calculations, we show the origin of these conclusions and essentially derive 
them by developing and using semi-analytical approximations. The expected amplitudes, shapes, 
oscillatory features, etc. of all correlation functions as functions of ^ are under analytical control. 
The central point of this analysis is the TE correlation function. We show that, at lower multipoles 
^, the TE correlation function must be negative (anticorrelation), if it is induced by gravitational 
waves, and positive, if it is induced by density perturbations. We argue that this difference in sign 
of TE correlations can be a signature more valuable observationally than the presence or absence 
of the BB correlations. This is because the TE signal is about two orders of magnitude larger 
than the expected BB signal and is much easier to measure. We summarize the competing effects 
of density perturbations in Appendix IdI 

Theoretical findings are compared with observations in Sec. IVIIII In the context of relic gravi- 
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tational waves it is especially important that the WMAP team |9| stresses (even if for a different 
reason) the actual detection of the TE anticorrelation near £ ~ 30. We show that this is possible 
only in the presence of a significant amount of relic gravitational waves (within the framework of 
all other common assumptions). We analyze the CMB data and suggest models with significant 
amounts of gravitational waves that are consistent with the measured TT, TE and EE correlation 
functions. Our final conclusion is that there is evidence of the presence of relic gravitational waves 
in the already available CMB data, and further study of the TE correlation at lower ts has the 
potential of a firm positive answer. 

II. GENERAL PROPERTIES OF RELIC GRAVITATIONAL WAVES 

Here we summarize some properties of cosmological perturbations of quantum-mechanical ori- 
gin. We will need formulas from this summary for our further calculations. 

A. Basic definitions 

As usual (for more details, see, for example, Q|), we write the perturbed gravitational field of 
a flat FLRW universe in the form 

ds'^ = -c^dt^ + a^{t){6ij + hij)dx'dx^ = o?{'q) [-d-rf + {6ij + hij)dx'dx^] . (1) 

We denote the present moment of time hy rj = r]^ and define it by the observed quantities, for 
example, by today's value of the Hubble parameter Hq = H{r]ji). In addition, we take the present 
day value of the scale factor to be a(r?ij) = 21h, where Ih = c/Hq. 

The functions hij (7/,x) are expanded over spatial Fourier harmonics e^*"'^, where n is a dimen- 
sionless time independent wave-vector. The wavenumber n is n = {dijn'n^Y^'^ . The wavelength 
A, measured in units of laboratory standards, is related to n by A = 2-Ka/n. The waves whose 
wavelengths today are equal to today's Hubble radius carry the wavenumber tih = 47r. Shorter 
waves have larger n's and longer waves have smaller n's. 

The often used dimensional wavenumber defined by /c = 2ti /\[r]{i) in terms of today's wave- 
length X{r]ji), is related to n by a simple formula 

k = ^^n (1.66 X 10"'' h) Mpc"^ 

The expansion of the field hij {rj, x) over Fourier components n requires a specification of polar- 

s 

ization tensors Pij (n) (s = 1, 2). They have different forms depending on whether the functions 
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hij(rj,x) represent gravitational waves, rotational perturbations, or density perturbations. 

In the case of gravitational waves, two independent linear polarization states can be described 
by two real polarization tensors 

1 2 

Pij (n) = klj - rmmj, Pij (n) = kmj + rmlj, (2) 

where spatial vectors (1, m, n/n) are unit and mutually orthogonal vectors. The polarization tensors 
((^j) satisfy the conditions 

Pij5'^ = 0, kjn' = 0, iijp'^ = 25s's- (3) 

1 2 

The eigenvectors of Pij are li and mj, whereas the eigenvectors of Pij are li + nii and U — rrii. In 
both cases, the first eigenvector has the eigenvalue +1, whereas the second eigenvector has the 
eigenvalue —1. 

In terms of spherical coordinates 0, (/), we choose for (1, m, n/n) the right-handed triplet: 

1 = (cos cos cos6'sin</), — sin0), 

m=(— sin0, COS0, 0), (4) 

n/n = (sin cos 0, sin sin 0, cosO). 
The vector 1 points along a meridian in the direction of increasing 0, while the vector m points 
along a parallel in the direction of increasing (j). With this specification, polarization tensors @ 
will be called the '+' and 'x' polarisations. The eigenvectors of '+' polarization correspond to 
north-south and east- west directions, whereas the ' x ' polarization describes the directions rotated 
by 45°. 

With a fixed n, the choice of vectors 1, m given by Eq. (0]) is not unique. The vectors can be 
subject to continuous and discrete transformations. The continuous transformation is performed 
by a rotation of the pair 1, m in the plane orthogonal to n: 

1' = Icos^ -|- msin-;/', m' = — Isin^/; -|- mcosV', (5) 

where V' is an arbitrary angle. The discrete transformation is described by the flips of the 1, m 
vectors: 

1' = —1, m' = m or l' = 1, m' = — m. (6) 
When (jSJ is applied, polarization tensors transform as 
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1 2 



p',j (n) = I'il'j - m'im'j =P,j (n) cos 2^+ P,j (n) sin 2^^, 

2 12 ^ ' 

p'ij (n) = I'-m'j + m[l'- = - Pij (n) sin 2^--^ Pij (n) cos 2ip, 



and when (jEJ is applied they transform as 

1 1 2 2 

P'y (n) =Pi, (n), p'ij (n) = - Pij (n). (8) 

Later in this section and in Appendix ^ we are discussing the conditions under which the 
averaged observational properties of a random field hij{n, rf) are symmetric with respect to rotations 
around the axis n/n and with respect to mirror reflections of the axes. Formally, this is expressed as 
the requirement of symmetry of the g.w. field correlation functions with respect to transformations 
© and ©. 

In this paper we will also be dealing with density perturbations. In this case, the polarization 
tensors are 

These polarization tensors satisfy the last of the conditions ©■ 

In the rigorous quantum-mechanical version of the theory, the functions hij are quantum- 
mechanical operators. We write them in the form: 



+00 



kj (n) hn iv) e'"''' Cn +Pij (n)L (f?)e ^"''^ 



(10) 



where the annihilation and creation operators, Cn and c^, satisfy the relationships 

[c„, v]= W('^(n-n'), c„|0)=0. (11) 

The initial vacuum state |0) of perturbations is defined at some moment of time rjQ in the remote 
past, long before the onset of the process of superadiabatic amplification. This quantum state 
is maintained (in the Heisenberg picture) until now. For gravitational waves, the normalization 



constant is C = VlGvr/p;. 

The relationships 1)11(1 determine the expectation values and correlation functions of cosmolog- 
ical perturbations themselves, and also of the CMB's temperature and polarization anisotropies 
caused by these cosmological perturbations. In particular, the variance of metric perturbations is 
given by 



oo 

-'2 



f 2 r s dn 

(0| /..,(r/,x)/.*^(r/,x) |0) = ^ / | hn (r/)!^-. (12) 

i s=l,2 

The quantity 

h\n,r^) = f^n' ^ | L (^)P = ^ E I ^ (l^) 

^ s=l,2 s=l,2 



is called the metric power spectrum. Note that we have introduced 

h{n,ri) = -nhn{ri). (14) 
vr 

The quantity (|1H|) gives the mean-square value of the gravitational field perturbations in a logarith- 
mic interval of n. The spectrum of the root-mean-square (rms) amplitude /i(n, rj) is determined by 

___ s 

the square root of Eq. (|13|) . Having evolved the classical mode functions hn {'>]) up to some arbi- 
trary instant of time ij (for instance, today rj = rj^i) one can find the power spectrum h{n, rf) at that 
instant of time. For the today's spectrum in terms of frequency u measured in Hz, v = uHq/Att, 
we use the notation hymsi^)- 

In our further applications we will also need the power spectrum of the first time-derivative of 
metric perturbations: 



^ ' drj drj I ^ 2 y ^f^^ 



oo 

dn 







d h (n, T]) 



drj 



(15) 

n 



To simplify calculations, in the rest of the paper we will be using a 'classical' version of the 
theory, whereby the quantum-mechanical operators and are treated as classical random 
complex numbers Cn and c^. It is assumed that they satisfy the relationships analogous to Hll() : 

(Cn) = (c^,) = 0, (SnC„,) = (S„c„,) = dss'S^'^Hn - n'), {Lcn') = iKK') = 0, (16) 

where the averaging is performed over the ensemble of all possible realizations of the random field 

m- 

The relationships 1)16(1 are the only statistical assumptions that we make. They fully determine 
all the expectations values and correlation functions that we will calculate, both for cosmological 
perturbations and for the induced CMB fluctuations. For example, the metric power spectrum 
(|13|) follows now from the calculation: 

oo 

1 {h,,{v,x)h^^{v,x)) = I L m'^. (17) 



s=l,2 





s 

The quantities | hn {vn ^-re responsible for the magnitude of the mean-square fluctuations 
of the field in the corresponding polarization states s. In general, the assumption of statistical 
independence of two linear polarization components in one polarization basis is not equivalent to 
this assumption in another polarization basis. As we show in Appendix statistical properties 
are independent of the basis (i.e. independent of tp, Eq. ©), if the condition 

I hn (ry)P = I hn (r?)!' (18) 



is satisfied. As for the discrete transformations Q, they leave the g.w. field correlation functions 
unchanged. 

In our further discussion of the CMB polarization it will be convenient to use also the expansion 
of hij over circular, rather than linear, polarization states. In terms of definitions ((21), the left 
(s = L) and right (s = R) circular polarization states are described by the complex polarization 
tensors 

L * R R * L ^ ' 

Pij Piji Pij Piji 

satisfying the conditions (for s = L,R) 

kjS'^=0, Pijn' = 0, iijp'^* = 26s's. (20) 



A continuous transformation ((SJ brings the tensors H19|) to the form 



^ 4., e-^'^ P\,=i,e^'t (21) 



Functions transforming accor ding to the rule ((21]) are called the spin-weighted functions of spin +2 
and spin -2, respectively A discrete transformation ^ applied to (fT^ interchanges 

the left and right polarisation states: 

P'ij=Pij, P'ij=kj. (22) 

The assumption of statistical independence of two linear polarisation states is, in general, not 
equivalent to the assumption of statistical independence of two circular polarization states. More- 
over, it is shown in Appendix^that symmetry between left and right is violated, unless 

I hn (r?)!' = I hn (r?)|^ (23) 

However, if conditions (|18|) . (|23|) are satisfied, statistical properties of the random g.w. field remain 
unchanged under transformations from one basis to another, including the transitions between 
linear and circular polarizations. The summation over s in the power spectra such as ()17|) can be 

s 

replaced by the multiplicative factor 2. Moreover, the mode functions hn {f]) ^'^^ t'^o independent 
polarization states will be equal up to a constant complex factor e*°. This factor can be incor- 
porated in the redefinition of random coefficients without violating the statistical assumptions 
(|16j) . After this, the index s over the mode functions can be dropped: 

L (r?) = hniv)- (24) 



There is no special reason for the quantum-mechanical generating mechanism to prefer one 
polarization state over another. It is natural to assume that the conditions (|18|1 . (|23|) hold true for 
relic gravitational waves, but in general they could be violated. In calculations below we often use 
these equalities, but we do not enforce them without warning. 



B. Mode functions and power spectra 

s 

The perturbed Einstein equations give rise to the g.w. equation for the mode functions hn (v)- 
This equation can be transformed to the equation for a parametrically disturbed oscillator ^J: 



An + A„ 



a" 



0, (25) 



where An {f]) = ^(^) hn (??) and ' = d/d-q = {a/c)d/dt. (In this paper, we ignore anisotropic 
stresses. For the most recent account of this subject, which includes earlier references, see 
Clearly, the behaviour of the mode functions depend on the gravitational 'pumping' field a(ry), 
regardless of the physical nature of the matter sources driving the cosmological scale factor a(r/). 
Observational data about relic gravitational waves allow us to make direct inferences about H{rj) 
and a(?7) and it is only through extra assumptions that we can make inferences about such 
things as, say, the scalar field potential (if it is relevant at all). 



Previous analytical calculations (see, for example. 
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and references there) were based on 



models where a{'q) consists of pieces of power- law evolution 

«(r/)=^o|r/|i+^ (26) 

s f> 

where lo and j3 are constants. The functions a(ry), a'(r/), hn (??)i h'n (??) were continuously joined 
at the transition points between various power-law eras. 

It is often claimed in the literature that this method of joining the solutions is unreliable, unless 
the wavelength is "much longer than the time taken for the transition to take place" . Specifically, 
it is claimed that the joining procedure leads to huge errors in the g.w. power spectrum for short 
waves. It is important to show that these claims are incorrect. As an example, we will consider 
the transition between the radiation-dominated and matter-dominated eras. 

The exact scale factor, which accounts for the simultaneous presence of matter, oc a(r/)~^, 
and radiation, pj oc a(ry)~^, has the form 



ai,) = 2ln{\±^] V h + , (27) 



2 ~\~ Z(,q J \ 1 -|- Z(^q 



10 



where z^q is the redshift of the era of equahty of energy densities in matter and radiation Pm{zeq) 
For this model, the values of parameter tj at equality and today are given by the expressions 



^eq 



1 -\- Z^q 1 -j- Z^q 1 ~1~ Z^q 

The current observations favour the value Zeq^3x 10^ Q- 

The piece-wise approximation to the scale factor (|27|) has the form 

aiv) = r], r] < rjeq] a{r]) = 21h {v + Veqf , r/ > Veq, (28) 

Y -L + Zgq 



where rj^q and rm for this joined a(r/) are given by r^eg = l/2y/l + Z(,q and = 1 — 1/2-^/1 + Zeg. 
It can be seen from ()27() and H28() that the relative difference between the two scale factors is very 
small in the deep radiation-dominated era, r/ ^ l/^Ze<j) and in the deep matter-dominated era, 
r/ ^ 1/ ^Zeq- But the difference reaches about 25% at times near equality. 

The initial conditions for the g.w. equation (|25() are the same in the models ()27p . ()28() and 
they are determined by quantum-mechanical assumptions at the stage (which we call the i-stage) 
preceding the radiation epoch. The i-stage has finished, and the radiation-dominated stage has 
started, at some rji with a redshift Zi. The numerical value of Zi should be somewhere near 10^^ 
(see below). 

The initial conditions at the radiation-dominated stage can be specified at that early time r/j or, 
in practice, for numerical calculations, at much latter time, as long as the appropriate g.w. solution 
is taken as 

s s 

Pn {ff) = —2iB sinnr], p^ (?]) = —2inB cos nr], (29) 

where 



B = F{p) 



+ Zeq 
l + Zi 



and F[[3) is a slowly varying function of the constant parameter (3, |i<'(— 2)| = 2. Parameter (5 
describes the power-law evolution at the i-stage and determines the primordial spectral index n, 
n = 2/3 + 5. In particular, (3 = —2 corresponds to the flat (scale-invariant) primordial spectrum 
n = 1. 
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For numerical calculations we use the constant B in the form: 



B = 2 



47ry/l + z, 



1 + Zi 



■eq 



) 




(30) 



The wave-equation (|25() with the scale factor (|27() can not be solved in elementary functions. 
However, it can be solved numerically using the initial data (|29|1 . Concretely, we have imposed 
the initial data H29|) at r]r = 10^^, which corresponds to the redshift z ~ 3 x lO'^, and have chosen 
Zeq = 6 X 10^ for illustration. Numerical solutions for hn{ri) / hn{rjr) characterized by different 
wavenumbers n are shown by solid lines in Fig. ^ We should compare these solutions with the 
joined solutions found on the joined evolution ()28j) for the same wavenumbers. 

The piece-wise scale factor (|28() allows one to write down the piece-wise analytical solutions to 
the g.w. equation (^5]) : 



where 7/2 = nrjeq. In Fig. ^we show the joined solutions ()31() by dotted curves. 

It is clear from Fig. ^ that the g.w. solutions as functions of r] are pretty much similar in the 
two models. The solid and dotted curves are slightly different near equality (where the relative 
difference between the scale factors is noticeably large) and only for modes that entered the Hubble 
radius around or before equality. Moreover, the g.w. amplitudes of the modes that entered the 
Hubble radius before equality gradually equalize in course of later evolution. There is nothing like 
a huge underestimation or overestimation of the high-frequency g.w. power that was alleged to 
happen in the joined model. 

It is easy to understand these features. Let us start from waves that entered the Hubble radius 
well after equality, i.e. waves with wavenumbers n/riH ^ y^l + z^g. As long as these waves are 
outside the Hubble radius, their amplitudes remain constant and equal in the two models, despite 
the fact that the scale factors are somewhat different near equality. The waves start oscillating in 
the regime where the relative difference between (|27() and H28p is small, and therefore the mode 
functions, as functions of r], coincide. 




^Jn{r] + 'qeq) [A„ J3/2 ( n{ri + rjeq)) - iBnJ-3/2 ( ri{r] + 7?eg))] , V > Veq 



(31) 
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1.2 




time ri 

FIG. 1: The g.w. mode functions hn{vi) / hn{rir) in a matter-radiation universe. The solid curves are 
numerically calculated solutions on the scale factor (|27|l . while the dotted curves are analytical solutions on 
the scale factor (I28|l . 

The waves with wavenumbers njnu ^ ^\ + Zgg enter the Hubble radius well before the equal- 
ity. They oscillate in the WKB regime according to the law /i„(r/) oc e"^^"^ / a(r]) , having started 
with equal amplitudes in the two models. Near equality, the mode functions are different in the 
two models, as much as the scale factors are different. But the relative difference between ()27() and 
H28|) decreases with time, and therefore the mode functions in the two models gradually equalize. 
The amplitudes of these mode functions, as well as the scale factors H28() . are exactly equal 

today. The only difference between these high-frequency mode functions is in phase, that is, in 
different numbers of cycles that they experienced by today. This is because the moment of time 
defined as 'today' in the two models is given, in terms of the common parameter r], by slightly 
different values of r]ji. 

Finally, for intermediate wavenumbers n/nn ~ y^l + Zeq, the modes enter the Hubble radius 
when the scale factors differ the most. Therefore, they start oscillating with somewhat different 
amplitudes. The difference between these mode functions is noticeable by the redshift of decoupling 
Zdec (characterized by somewhat different values of rjdec), as shown in Fig.^ The difference survives 
until today, making the graph for the g.w. power spectrum (|1.S|) a little smoother (in comparison 
with that derived from the joined model) in the region of frequencies 10~^^ Hz that correspond to 
the era of equality. 
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frequency 



FIG. 2: The present-day spectra for hrms{v) and Vlgyj{v). The soHd hnes correspond to the primordial 
spectral index (3 = —1.9, i.e. n = 1.2, while the dashed lines are for j3 ~ —2, i.e. n = 1. 

Having demonstrated that the use of joined analytical solutions is well justified, and the previ- 
ously plotted graphs for hrmsi'^) ^-nd r2g^(i/) are essentially correct, we shall now exhibit the more 
accurate graphs based on numerical calculations with the initial conditions ()29() . H3U|) . We adopt 
Zeg = 3 X 10^ and Hq = 75^ /Mpc We adjust the remaining free constant Zi in such a way 
that the temperature correlation function £{i + l)Ci/2'iT at ^ = 2 is equal to 211 /zK^ This 
requires us to take Zi = 1.0 x 10^^ for f3 = —2 and Zi = 2.4 x 10^° for /3 = —1.9. 

The graphs for today's spectra, normalised as stated above, are shown in Fig. [2 In order to 
keep the figure uncluttered, the oscillations of hrms{t^) are shown only at low frequencies, while the 
oscillations of Qg^^i') are not shown at all. We have also indicated possible detection techniques 
in various frequency bands. 

The function Qgw{u) is the spectral value of the cosmological parameter [3]: 

\lgw{l'l,l'2) = — = — / Pgw[v) = / ilgw{v) , 

Pc Pc J Ul ^ Jui ^ 

that is, 

^gw{l^) = —Pgw{l^)- 
Pc 

Using the high-frequency definition of Pgw{^) (valid only for waves which are significantly shorter 
than Ih) we derive 



^gM = ^^hl^,[u)(^\\ (32) 



3 K'^hJ 
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It is this definition of 57^^(1/) that is used for drawing the curves in Fig. [21 

We have to warn the reader that a great deal of hterature on stochastic g.w. backgrounds uses 
the incorrect definition 



1 dpgu,{l>) 



Pc d Inu ' 

which suggests that the 17^^(1/) parameter is zero if the g.w. spectral energy density Pgw{^') is 
frequency-independent, regardless of the numerical value of Pgwi^^)- Then, from this incorrect 
definition, a formula similar to Eq. (|32() is often being derived by making further compensating 
errors. 

The higher- frequency part of hrms{v) is relevant for direct searches for relic gravitational waves, 
while the lower-frequency part is relevant to the CMB calculations that we turn to in the next 
section. The direct and indirect methods of detecting relic gravitational waves are considered in 



a large number of papers (see, for example, references 
though we disagree with some of them). 
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26 



27 



even 



It is important to keep in mind that according to the infiationary theory the dashed lines for 
^gwi^^) and hrms{i^) in Fig. I^lshould be at a zero level, because they describe the g.w. background 
with a flat primordial spectrum /? = —2, n = 1. The 'consistency relation' of the inflationary 
theory demands r = 0, i.e. vanishingly small g.w. background, in the limit n = 1. 



III. EQUATIONS OF RADIATIVE TRANSFER 
A. Characterization of a radiation field 



The radiation field is usually characterized by the four Stokes parameters (/, Q, U, V), 23. l3C | . 
The parameter / is the total intensity of radiation, Q and U describe the magnitude and direction 
of linear polarization, and V is the circular polarization. The Stokes parameters can be viewed 
as functions of photons' coordinates and momenta Since photons propagate with the 

speed of light c, the momenta satisfy the condition PaP"' = 0. Thus, the Stokes parameters are 
functions of (t, x*, u, e*), where v is the photon's frequency, and e* is a unit vector in the direction 
of observation (opposite to the photon's propagation). 

In a given space-time point (t, x^), the Stokes parameters are functions of i', 0, (p, where 6, (p are 
coordinates on a unit sphere: 

da^ = gahdx^dJ' = dO'^ + sin^ Odc/)'^. (33) 
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The radial direction is the direction of observation. 

The Stokes parameters are components of the polarization tensor Pat jsol which can be written 

as 

1 / I + Q - (U-iV)sm9 \ 

\-{U + iV)sme {I-Q)sm^ej 

(We do not indicate the dependence of Stokes parameters on ly.) 

The symmetry of Pab with respect to rotations around the direction of observation requires the 
vanishing of linear polarization Q = 0, U = 0, but the circular polarization V can be present. The 
symmetry of Pab with respect to coordinate reflections in the observation plane requires = 0, 
but the linear polarization can be present. The first symmetry means that the readings of a linear 
polarimeter are the same when it is rotated in the observation plane. The second symmetry means 
that the readings of the left-handed and right-handed circular polarimeters are the same. (Compare 
with the discussion on gravitational waves in Sec. Ill Al ) 

Under arbitrary transformations of 9,(j), the components of Pab{0,(t>) transform as components 
of a tensor, but some quantities remain invariant. We want to build linear invariants from Pab 
and its derivatives, using the metric tensor Qabid, (f) and a completely antisymmetric pseudo-tensor 
e'^^(0,0). 



^ab 



- sin"^ 
sin-i e 



The first two invariants are easy to build: 

Iie,^)=g''\9,cP)Pabi9A), V{e,^) = ie''\e,^)Pabi9A)- (35) 



Then, it is convenient to single out the trace and antisymmetric parts of Pab, and introduce the 

ab 

1 



symmetric trace- free (STF) part Pa^^ ' 



Q -UsmO 



tdSTF 
^ab 



-Usm9 -Qsm'^9 



Clearly, the construction of other linear invariants requires the use of covariant derivatives of 
the tensor Paif^- There is no invariants that can be built from the first derivatives Pa^^ , so we 
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need to go to the second derivatives. One can check that there are only two linearly independent 
invariants that can be built from the second derivatives: 



:b:d 



(36) 



The quantities / and E are scalar s, while V and B are pseudoscalars. V and B change sign under 
flips of directions (coordinate transformations with negative determinants). This is also seen from 
the fact that their construction involves the pseudo-tensor e^fe. 

The invariant quantities (/, E, B), as functions oi 9,(j), can be expanded over ordinary spher- 
ical harmonics Yii„i{0 ^ (j)) ^ = 



1(0 A 
v{e,(i) 

E{9,4> 
B{0,(l) 



oo i 

f=0 m=-l 

oo I 



t=Q m=-e 
oo i 




(37a) 
(37b) 
(37c) 
(37d) 



The set of multipole coefficients {aj^, a^^, af^, af^) completely characterizes the radiation field. 
We will use these quantities in our further discussion. 

To make contact with previous work, we note that the multipole coefficients CLf^,af^ can also 
be expressed in terms of the tensor Pab itself, rather than its derivatives. This is possible because 
one can interchange the order of differentiation under the integrals that define a^, af^ in terms of 
the right hand side (r.h.s.) of Eq. This leads to the appearance of the spin-weighted spherical 



harmonics or tensor spherical harmonics 
written as 



12 



3 



14 



32l |. For example, the tensor Pab can be 



^ oo i 



le 



V 

ab'^im 



£=0 m 

oo 



+ 



i=2 m=-l 



where Y,9 



{£m)afe(^' and the "gradient" and "curl" tensor spherical harmonics 

forming a set of orthonormal functions for STF tensors The invariants E and B can also be 

n 

written in terms of the spin raising and lowering operators 9 and 3 |31l |: 



1 



E = -- [§2 (Q + iu) + g2 (Q _ iu)] , B = -[n^{Q + iu) - g2 (g - iu)] . 
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The quantities E and B are called the E (or "gradient") and B (or "curl") modes of polarization. 
The ^-dependent numerical coefficients in H37c|) and ()37d|) were introduced in order to make the 



definitions of this paper fully consistent with the previous literature 



31 



32 



B. Radiative transfer in a perturbed universe 



We need to work out the radiative transfer equation in a slightly perturbed universe described 
by Eq. The Thompson scattering of initially unpolarized light cannot gene rate circular polar- 



29 



33 



34J . we shall write the 



n 



1^ 

2/11/3 



I-Q 
-2U 



(38) 



ization, so we shall not consider the V Stokes parameter. Following 
radiative transfer equation in terms of a 3-component quantity (symbolic vector) n(x",p"). The 
components (ni,n2,ri3) are related to the Stokes parameters by 

(l + n\ 

1 C2 

where h is the Planck constant. The quantities fii,h2, (fii + n2 + n-3)/2 are the numbers of photons 
of frequency z/ coming from the direction z and passing through a slit oriented, respectively, in the 
directions x, y, and the bisecting direction between x and y. 

The equation of radiative transfer can be treated as a Boltzmann equation in a phase space. 
The general form of this equation is as follows 



35| 



as 



(39) 



where s is a parameter along the world-line of a photon, ^ is a total derivative along this world-line, 
and C is a collision term. We shall explain each term of this equation separately. 
The total derivative in Eq. (|39() reads: 



Dh 

ds 



dx" d dp"' d 
+ 



n. 



ds dx" ds dp" 

where dx" /ds and dp" /ds are determined by the light-like geodesic world-line. 



dx" 
ds 



P , 



dp" 
ds 



r" rJ^n'y n ort"n 
A P ) QafiP P 



0. 



(40) 



(41) 



Strictly speaking, the square bracket in Eq. (|4n|) should also include the additive matrix term R. 
This term is responsible for the rotation of polarization axes that may take place in course of 



parallel transport along the photon's geodesic line 



3q |. In the perturbation theory that we are 
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working with, this matrix does not enter the equations in the zeroth and first order approximations 
[3^, and therefore we neglect R. 

In our problem, the collision term C describes the Thompson scattering of light on free (not 
combined in atoms) electrons 29(. We assume that the electrons are at rest with respect to one 
of synchronous coordinate systems (^. We work with this coordinate system, so that it is not 
only synchronous but also 'comoving' with the electrons. (This choice is always possible when 
the functions hij in ^ are gravitational waves. Certain complications in the case of density 
perturbations will be considered later, Appendix IdI) Thus, the collision term C is given by the 
expression 



C [n] = -arNeix'^) (^^^ h{t, x\ v,e,<P)-^ j dn'P{e, 0; 9', (l)')h{t, x\ v, e\ 0') 



(42) 



where ax = 6.65 • 10^^^ cm^ is the Thompson cross section, N^, is the density of free electrons, and 
F{6,(j);6' is the Chandrasekhar scattering matrix. (The explicit form of the scattering matrix 
is discussed in Appendix lO) The factor cdt/ds arises because of our use of the element ds, instead 
of cdt, in the left hand side (l.h.s.) of Eq. H39() . In accord with the meaning of the scattering term, 
the quantity aTNe{x°'){cdt / ds) is the averaged number of electrons that could participate in the 
scattering process when the photon traversed the element ds along its path. 

Let us now write down the equations of radiative transfer in the presence of the gravitational 
field (^. First, we write down the equations for the light-like geodesic line a;"(s) = (77(5), a;*(s)): 



r, dn u dx* v ,• du 

p^ = -r = —, P =-r = — e ' -r 

ds ca ds ca ds 



{5ij + hij) e'e^ = 1. 



1 da 1 A dhij 

1 — e e — 

adrj 2 di] 



We do not need the expression for de^ /ds, because it is a first-order (in terms of metric perturba- 
tions) quantity, and this quantity enters the equations of radiative transfer only in products with 
other first-order terms. We neglect such second-order corrections. 
Second, we write for the 'vector' n: 

n = n(o)+n«, (44) 

where 

n(0) is 

the zeroth order solution, and fi^^) is the first order correction arising because of 
the presence of metric perturbations. We shall now formulate the equation for n^^^, taking into 
account the zero-order solution to Eq. (jSH). 
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In the zero-order approximation we assume that hij = and that the radiation field is fully 
homogeneous, isotropic, and unpolarized. Therefore, 



n 



W =nW(r?,i/)u, 



(45) 



where 

u = 1 

vv 

Since the scattering matrix P does not couple to the radiation field if it has no quadrupole 
anisotropy, the collision term H42|) vanishes, C [n'-*'-*] = 0, and the equation for n^^\rj,i') reads 

9n(°) u da an(0) _ ^ 
dr] a dr] dv 

The general solution to this equation is n^'^'> = uq {1/0(7])), which makes it convenient to use a new 
variable 

D = va{rj). 

In the zero-order approximation, v = const /a{rj) and therefore i> is a constant along the light ray. 

We are now in a position to write down the first-order approximation to the Boltzmann equation 
((Snj- We take {r],x\i>,e^) as independent variables, i.e. n^^^ = no{h'), n.^^^ = n^-^^ (r/, x*, i>, e*), and 
use the identity 

dp" d di> d ^de^ d 
ds dp" ds du ds de^ 

in the first-order approximation to ()4U() . Taking also into account the geodesic equation H43|l we 

arrive at the equation 



h e : ve. e' 

dr] dx^ 2 df] dv 



^ = C 

ds 



fid) 



Introducing new notations g(ry) = aTOi{r])N e,(r]) and f(y) = (?ln no/(? In t' (the astrophysical 
meaning and numerical values of the functions qij]) and f(y) are discussed in Appendix El) we 
write down the final form of the transfer equation: 



d ^ ^^^^ _^ g« ^ 
dr] dx^ 
f{i')no{i>) j ,dhij , 

2 d^ 



u + g(r?)^ y P(e^e'^)n(i)(r?,x\i>,e'-'). 



(46) 
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It is seen from Eq. (|46() that the 'source' for the generation of n*-^-* consists of two terms on 
the r.h.s. of this equation. First, it is the gravitational field perturbation hij, participating in the 
combination e^e^dhij/drj. It directly generates a structure proportional to u, i.e. a variation in 
the / Stokes parameter and a temperature anisotropy. In this process, a quadrupole component of 
the temperature anisotropy necessarily arises due to the presence of the term e^d/dx^, even if the 
above-mentioned combination itself does not have angular dependence. The second term on the 
r.h.s. of Eq. H46|) generates polarization, i.e. a structure different from u. This happens because 
of the mixing of different components of h^^\ including those proportional to u, in the product 
term Pn^^^. In other words, polarization is generated by the scattering of anisotropic radiation 
field. Clearly, polarization is generated only in the intervals of time when q{r]) ^ 0, i.e. when free 
electrons are available for the Thompson scattering (see, for example, js^). 



C. The radiative transfer equations for a single gravitational wave 

We work with a random gravitational field hij expanded over spatial Fourier components (|1U() . 
It is convenient to make similar expansion for the quantities n^^\r],x^,i',e'^). Since Eq. (|46() is 
linear, the Fourier components of n^^^ inherit the same random coefficients Cn that enter Eq. ()1U() : 

+ CX) 



(27r)3/2 



/2n 



s=1.2 



(47) 



Equation (|46jl for a particular Fourier component takes the form: 



dr] 



+ Qiv) + ^eVj 



= l^^e^e^ P, („)%Mu + Mfdn' P(e^ e'^)n« (., ^, e'^- (48) 
2 dr] 47r J 

To simplify technical details, we start with a single gravitational wave propagating exactly in 
the direction of z in terms of definitions 0, @, (|19jl . The coordinate system, associated with 
the wave, is specified by = 0, = 0. This simplifies the polarization tensors H19|) and makes 
them constant matrices. At the same time, the observation direction is arbitrary and is defined 
by e* = (sin0cos(/), sin sine/;, cos 9). We consider circularly polarized states with s = 1 = L, 
s = 2 = R. Then, we find 



e'ei Pij (n) = (1 - ]i^)e 



(49) 
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where ^ = cos 9 and the it signs correspond to s = L,R, respectively. This simpUfication of the 
angular dependence is possible only for one Fourier component, but not for all of them together. 
We shall still need the results for a wave propagating in an arbitrary direction. The necessary 
generalization will be done in Sec. IV Bl 

The ziz2(j) angular dependence of the source term in Eq. (|48|) generates the ±2(1) angular depen- 
dence in the solution 23, js^. We show in Appendix [0 that the terms in h^]i{r],i>, with any 
other ^-dependence satisfy homogeneous differential equations, and therefore they vanish if they 
were not present initially (which we always assume). Similarly, the u dependence of the solution 
can be factored out. Finally, we show in Appendix O that one linear combination of the three 
components of fin ^ always satisfies a homogeneous equation and therefore vanishes at zero initial 
data. Thus, the problem of solving Eq. (|15|l reduces to the problem of finding two functions of the 
arguments r], fj,. 

Explicitly, we can now write 

(l\ ( (1 + /^') ^ 



1 



(50) 



Clearly, function a is responsible for temperature anisotropy (I Stokes parameter), while function 
f3 is responsible for polarization {Q and U Stokes parameters). 



Temporarily dropping out the labels n, s and introducing the auxi^ 



a{r],fi) + P{ri,ij), we get from Eqs. (|i8|) . (|50|) a pair of coupled equations 34l |: 



iary function £,{rj,fi) 



dr] 



where 



1 



dr] 



(l + ^'Y/3(7?,^')-2(l-^")'^(^'/^') 



(51) 
(52) 

(53) 



IV. RADIATIVE TRANSFER EQUATIONS AS A SINGLE INTEGRAL EQUATION 

In some previous studies [^, 0, 0|, equations (|^ . are being solved by first expanding 
the /^-dependence of /3 and ^ in terms of Legendre polynomials. This generates an infinite series of 
coupled ordinary differential equations. Then, the series is being truncated at some order. 
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We go by a different road. We demonstrate that the problem can be reduced to a single 
mathematically consistent integral equation. There are technical and interpretational advantages 
in this approach. The integral equation enables us to derive physically transparent analytical 
solutions and make reliable estimates of the generated polarization. The numerical implementation 
of the integral equation considerably saves time and allows simple control of accuracy. 



A. Derivation of the integral equation 



In order to show that the solutions of Eqs. (|51jl . (|52|) for a{r],fj.) and (3(ri, fi) are completely 



determined by a single integral equation, we first introduce new quantities 



4C 



lb 



(54) 



H{r]) = e 



Solutions to Eqs. (j^T|) . can be written as 



-T(n) dhji]) 
dt] 



(55) 



(56) 



(57) 



Expression H56|) is a formal solution to Eq. H51|) in the sense that P{rj,fj,) is expressed in terms of 
$(ry) which itself depends on (3{r],fi) (see (|5:-{|) and (jSU)- 

We now put (|56|) and (|57j) into Eq. (|53|) to get a new formulation for X(r/) : 



1 r, 

Z{rj) = e^^'') j j dfidr]' 

-1 

Using the kernels K±{rj — rj'), 



(58) 



K±irj-V') = I rfMli^') 



(59) 



Eq. (|5S|) can be rewritten as 



X(r]) = e^('') j dri 




K^(j] - ) - -K^{r] - r]')H{r]') 



(60) 
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Multiplying both sides of this equality by (3/16)g(ry)e '^^'^^ and recalling the definition ((511) we 
arrive at a closed form equation for ^{r]): 

■n 

HV) = ^q{ri) I drj'^rj')K+{r, - r,') + F{rj), (61) 


where F(r]) is the known gravitational-field term given by the metric perturbations, 

F{ri) = -^Qiv) I dri'H{rj')K^{r] - r?'), (62) 



The derived equation (|61() for ^{rj) is the integral equation of Voltairre type. As soon as <&(??) is 
found from this equation, we can find /3(ry, /i) from Eq. (|56|) . Then, Eqs. (|57|) and (|56|) completely 
determine all the components of n^^^ according to Eq. (j50() . 

Clearly, we are mainly interested in temperature and polarization anisotropics seen at the 
present time rj = rjj^. Introducing ( = n{r]R — rj) and restoring the indices n and s, we obtain the 
present-day values of a and j3: 

VR. 

an,s{lj) = an,s{VR, fJ-) = y [Hn,s{v) - ^n,si'n)^ c"*'"^, (63a) 



Vr 

PnAf^) = PnAm,f^) = I dr, $„,,(r/)e-^'^^. (63b) 



The integrals can safely be taken from 77 = as the optical depth r quickly becomes very large in 
the early Universe, and the 'source' functions Hn^siv) ^^^^1 <I>„^s(r/) quickly vanish there. We will 
work with expressions (|63|) in our further calculations. 



B. Analytical solution to the integral equation 

The integral equation (j6H) can be solved analytically in the form of a series expansion. Although 
our graphs and physical conclusions in this paper are based on the exact numerical solution to Eq. 
(|61|) . it is important to have a simple analytical approximation to the exact numerical solution. 
We will show below why the infinite series can be accurately approximated by its first term and 
how this simplification helps in physical understanding of the derived numerical results. 

We start with the transformation of kernels ()59() of the integral equation ()61|) . Using the identity 
Ij^k^ix/i _ (^(}_li(lx)^e^^^ , the kernels can be written as 

+1 2 
KAV-V) = jd,{l±,'fe-^^^-^') =2(lT^) 
-1 
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where x = n{r] — rj'). Now, taking into account the expansion 



smx 



X 



E 

m=0 



7X , 



(2m + 1)!' 



the r.h.s. of Eq. (jblj) can be presented in the form of a series: 



^r^) = -q{r]) n^" / dr]' {rj - r^'f'' A+(m)$(r/') - \.{m)H{r]' 



where 



A+(m) 



(-ir 



(2m + 1)! 



m=0 



1-4 



(64) 



(m + 2) 



(2m + 3)(2m + 5) 



A-_(m) 



-1) 



(2m + l)!(2m + 3)(2m + 5) ' 



Since the r.h.s. of Eq. (|64|) is a series in even powers of the wavenumber n, the l.h.s. of the 
same equation can also be expanded in powers of v?'^: 



,2m 



(65) 



m=0 



Using expansion H65|) in both sides of Eq. H64|) we transform this equation to 
J]$(-)(ry)n2™= -g(r/) 



m=0 



^ ^ A+(m)n 

m=0 j=0 



2{m+j) 



j dri'^^^\ri'){'q-ri'] 


^ A_ (m)n2'" y dr]'H{ri'){r] - rjf' 



2m 



(66) 



The left side and the second term in the right side of Eq. (|66() are series in n?^, but the first sum 
on the r.h.s. of Eq. (|66|) is still a mixture of different powers. This sum can be rearranged to be 
manifestly a series in n^"*: 

V 



^ ^ A+(m)n2('^+j) / d7?'$(j)(r?')(r/ - r/'; 



2m 



m=0 j=0 



m=0 



2m 



f;A+(A:) '=)(r?')(7?-V; 
fc=o ;( 



2A: 



(67) 



According to Eq. ()66() we have to make equal the coefficients of terms with the same power n^™' 
in both sides of the equation. This produces a set of integral equations 



(68) 
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where functions S^'^\rj) depend only on the known function H{r]) and functions ^^"^ ''^v) Pi'^- 
sumed to be found from equations of previous orders: 

S(-"-^)(r,) = -^A_(m) j dT]' H{r]'){r]-i]'f"' 


V 

+^f;^+(^) / ^V' ^^"'-'Hv'Xv-vT- (69) 
k=i i 

The important advantage of the performed expansion in powers of n^™ is that Eq. (|68|) of any 
order m is now a self-contained analytically solvable integral equation. 

Exact solution to the integral equation H68|) is given by the formula 

v 

J df]' 


We can further simplify this formula. Taking an //-derivative of expression (|69|1 we find 



J\2m-1 



^ = -Y^H{rj)dom - uj-{m) J drj H{t] ){v-V J 





k=l 

where uj±{m) = 3mX±{k). Substituting this expression into Eq. H70|). we arrive at the final result 

V 







(71) 



+ 5]u;+(A;)cI>(— ^■)(^')^(fe)(^,r/' 

k=l 

where 

*W(^,^0 = / dv" e-i^-^^"'^'Hv" - VT'-'- (72) 
n' 

Functions ^ (^j^^[r],'q') depend only on the ionization history of the background cosmological 
model described by q{ri). These functions can be computed in advance. Complete determination 
of the function <I>(r?), Eq. requires only one integration by r/ at each level m in Eq. (|71|) . 

starting from m = 0. The zero-order term (^^^\'q) does not depend on functions '^(k) and is 
determined exclusively by H{i]), Eq. ()55p. The zero-order term can be presented as 

V 

(r?) = -^^9{V) I drj'^!^e-^o^(^'^'\ (73) 
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It is crucial to remember that the function $(ry), Eq. (|54() . always contains the narrow visibility 
function g(r]) (see Appendix IH)) . In particular, function ^^^\ri) is nonzero only for rj within 
the width of g{r]), and is proportional to this width. In the era of decoupling, we denote the 
characteristic width of g{r]) by Ary^ec- With Aijdec we associate the characteristic wavenumber n*: 

2tt 
= -i: • 

Numerically, Aj^^gc ~ 3 x 10~^ and n,,, ~ 2 x 10^. In what follows, we will be interested in CMB 
multipoles i < 10'^. They are mostly generated by perturbations with wavenumbers n < 10^. 
Therefore, for wavenumbers of interest, we regard n/n^: as a small parameter. 

We shall now show that <I>(*^)(r/) is the dominant term of the series ()65() . The next term, 
<I>(^)(?7)n^, is at least a factor (n/n^,)^ smaller than <I>(*^)(r/), and so on. The explicit expression for 
^^^^ is as follows 



Effectively, function (^^^\ri), in comparison with ^>W(7?), contains an extra factor ^'(^^^(r], r]'). Tak- 
ing into account the fact that the functions g{ri) and e localized in the interval of 
arguments not larger than Arjdec, this factor evaluates to a number not larger than (Ar/f^ec)^- 
Therefore, the m = 1 term in Eq. H65|) is at least a factor (n/n^)'^ smaller than the m = term. 

These analytical evaluations are confirmed by numerical analysis as shown in Fig. |5J The solid 
line shows the exact numerical solution found from Eq. ()61() for hn{ri) and q{rj) described in Sec. Ill Bl 
and Appendix^ respectively. The dashed line is plotted according to formula ((75|) . with the same 
/i„(r/) and q{r]). It is seen from Fig. EJthat the zero-order term <I>('^)(r/) is a good approximation. 
The deviations are significant, they reach (20-25)%, only for the largest wavenumbers n in the 
domain of our interest. 



■(1 



dr]' 



(l)e^(V)$(0)(^/) 



■^-('''^')*(i)(r?,r?' 



V. MULTIPOLE EXPANSION AND POWER SPECTRA OF THE RADIATION FIELD 
A. Multipole coefficients 

Having found ^n,si'n) for a single gravitational wave specified by Eq. H49|l one can find a and 
P functions according to Eqs. (jlvEj) . (|63b|) . Then, using Eqs. (jMl), (El), ^ and (jHEI), one 

can find the multipole coefficients {X = I, E, B) participating in the decompositions pT]) . 
Although this route has been traversed before |3l|, we have made independent calculations in a 
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FIG. 3: Function ^n{v) for different values oin. The solid line is exact numerical solution to Eq. The 
dashed line is the zero order approximation (|73|l . The dotted line shows the approximation (|91(l (see below). 
The g.w. mode functions are normalized such that hn{rir) = 1. 

more general arrangement. Formulas derived in this subsection are effectively a confirmation of 
the correctness of calculations in Ref. ,31] • 

First, we integrate over photon frequencies (fot the definition of 7 see Eq. dHll)) and arrive at 
the following expressions 



En,s (/", </>) 
Bn,s (/", <t>) 



-7 



(1 - f^') ((1 + /^') 1^ + 8^1^ + 12) f3n,s e 



(74a) 
(74b) 

(74c) 



where, as before, the upper and lower signs correspond to s = 1 = L and s = 2 = R, respectively. 
The ziz2(j) dependence in Eqs. (|74j) implies that only the m = ±2 multipoles are nonzero. 

Then, we integrate Eqs. (fTi]) over angular variables in order to find according to Eq. pTj) . 
Using the notations T£{(), Ei((), Bi{() for the functions arising in course of calculations (and 
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called multipole projection functions) 

m) = ^ 

EdO = 
BeiC) = 



' + 2)! MC) 
'-2)\ C ' 
l{l - 1 



c 



MO+je-iiO 



and replacing a(/x) and by their expressions (|63|) . we finally arrive at 

aj^{n,s) = i-iY~'^ (^52,mSi,s + S-.2,mS2,s^ af {n, s), 

af^{n,s) = {-iY~^ (^S2,mSi,s - 6-2,mh,s^ af {n, s), 



where 



Ti 

at 



{n,s) = 7 V47r(2£ + 1) / dr, (Hr^ - ^-n.s (??)) T>(C), 





VR 



\n,s) = ^^^T^(2i + 1) / $„,,(r/)^,(C), 





Vr 



\n,s) = 7 V4vr(2£ + 1) ^ dr] cI>„,,(r/)i?^(C). 



(75a) 
(75b) 
(75c) 

(76a) 
(76b) 
(76c) 

(77a) 
(77b) 
(77c) 



B. Superposition of gravitational waves with arbitrary wavevectors 

It is important to remember that the result (|76|) is valid only for a special wave, with the 
wavevector n oriented exactly along the coordinate axis z. Since the perturbed gravitational 
field is a random collection of waves with all possible wavevectors n, and we are interested in 
their summarized effect as seen in some fixed observational direction ^, </), we have to find the 
generalization of Eq. H76|) to an arbitrary wave, and then to sum them up. 

To find the effect of an arbitrary wave, there is no need to do new calculations. It is convenient 
to treat calculations in Sec. IV Al as done in a (primed) coordinate system specially adjusted to a 
given wave in such a manner that the wave propagates along n' = (0,0, n). The observational 
direction e* is characterized by 0',(f)'. The quantities X = I,E,B calculated in Sec. IV Al are 
functions of 9',(l)' expanded over Yim{d' , cp'), 



(78) 



=0 m= 
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where af^ is a set of coefficients (compare with Eq. p7])l: a^, [(£ + 2)!/(£ — 2)!]^''^ a^, 
[(£ + 2)!/(^-2)!]i/2af„. 

Now, imagine that this special (primed) coordinate system is rotated with respect to the ob- 
server's (unprimed) coordinate system by some Euler angles 



a = (t>n, P = 0n, 7 = 



(see, for example 



4l|). The same observational direction e* is now characterized by 9,(1), and the 



same wavevector n' is now characterized by the unit vector 

n = n/n = (sin 9^ cos sin^n sin (p^, cos ^n)- 

Obviously, the already calculated numerical values of the invariants X{9',(f)') do not depend on 
the rotation of the coordinate system. Being expressed in terms oi 9,(j), the invariants describe the 
effect produced by a wave with a given (arbitrary) unit wavevector n, as seen in the direction 9, cj). 

The transformation between coordinate systems {9' , cf)') and {9, (p) is accompanied by the trans- 
formation of spherical harmonics, 

e 

Yt^(9',4l)= Di,^^{h)Y,rn'{9,<t>), 

m'=-l 

where 

^m',m (n) = r)ln<,m {<Pn,9n,0) 

are the Wigner symbols [4^. Later, we will need their orthogonality relationship 

J dnDl^pih)Dil^in) = J^5u'Smm' , (79) 

where = sin 9nd9nd<pn- 

We can now rewrite Eq. ((TSj) in terms of 9, (p, and thus find the contribution of a single arbitrary 
Fourier component. 



00 e / e \ 
e=o m=-e \m'=-e J 



(80) 



The superposition of all Fourier components of the perturbed gravitational field gives, at the 
observer's position x = and at t] = rjji (see (|47p ) . the final result: 



,~oo 

X{9,<P) = / ^ E [XnAO,cP) +X*^A9,<P) t] . (81) 



-00 



s=l,2 
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From this expression, combined with Eq. (|8U() . one can read off the random multipole coefficients 
afj^ that participate in the expansions (|37|): 




(82) 

Since X is a real field, the multipole coefficients obey the reality conditions 

as is seen directly from H82() . The gravitational- wave nature of metric perturbations is encoded in 
concrete values of the coefficients (f77)) . (|^ . But in all other aspects the argumentation presented 
here is general. 



C. Angular power spectra for temperature and polarization anisotropies 

It follows from Eq. (|82|) that the statistical properties of the multipole coefficients are fully 
determined by the statistical properties of the gravitational field perturbations represented by the 
random coefficients c^. A particular realization of is responsible for the particular realization 
of af^ actually observed in the sky. Having derived the distribution function for from some 
fundamental considerations (for example, from the assumption of the initial quantum-mechanical 
vacuum state of perturbations) we could estimate the probability of the observed set af^ within the 
ensemble of all possible sets. We could also evaluate the inevitable uncertainty in the observational 
determination of the parameters of the underlying random process. This uncertainty is associated 
with the inherent absence of ergodicity of any random process on a 2-sphere (i.e. sky) I4J]. In this 
paper, however, we adopt a minimalistic approach; we postulate only the relationships H16|) and 
calculate only the quadratic correlation functions for a^. 

Clearly, the mean values of the multipole coefficients are zeros, 

To calculate the variances and cross-correlation functions, we have to form the products af^a^^, 
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and then take their statistical averages. First, we find 

f n^dndQ 



/ X* X' \ ^ f n dndVt y—^ ^—^ ^—^ 

(a^™«^'m') = 7^ y 2n 2^ 2^ 2^ 



s=\,2mY=-lra\=-l 

'4m^ {n, s)a^^, (n, s)D^^^^ (n)I)^/„/ (n) 

{ri, s)af^,^ (n, s)D'_^^„,^ {^)D% ^„,.^ (fi)] . (84) 

We now take into account the fact (compare with Eq. ()76p) that 

«ftni iji, s)a^^, (n, s) oc 5^^^/ . 

This property allows us to get rid of summation over m'l in Eq. H84j) . Then, we perform integration 
over dJl and use the orthogonality relationships (|79() . We finally arrive at 

{af^a^^,) = See' 5mm', (85) 

where 



Cf^' = 4^2(2^ + 1) J '^dnY^ [afm{n,s)afm{n,s) + afm{n,s)afm{n,s)j . (86) 



Other quadratic averages, such as {O'fm^^m')^ i^fm'^em')^ follow from and the reality condition 
(jHl. 

The angular correlation and cross-correlation functions of the fields /, E, B are directly express- 
ible in terms of Eq. (|8H|) . For example, 

(/(^i,<Ai)/(^2,<A2)) = r(<5) = ^_^cr^'£(cos5), 

£=0 

where 5 is the angular separation between the directions (0i,(/)i) and {02,4>2) on the sky. If the 
actually measured values of af^ represent a particular realization of a gaussian random process, 
quantities Cf^', constructed from the measured af^ according to the r.h.s. of Eq. (|HH|l. are the 
best unbiased estimates 

One can note that the final result (|85() , ()86() is the integral of individual contributions (|76() from 
single gravitational waves given in a special frame discussed in Sec. IV A[ However, one cannot 
jump directly from (|76|1 to (|86|) (which is often done in the literature). In general, Eq. (|86|) does 
not follow from Eq. ()76|). By calculations given in this subsection we have rigorously shown that 
Eq. (|H?t|) is justified only if special statistical assumptions (|TH|) are adopted. 

One can also note that the correlation functions containing the label B once, i.e. CP and Cf ^, 
vanish if the extra assumptions ((TH|l. are made. Indeed, under these assumptions one can use 
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FIG. 4: The 'source' functions Hn.s{f]) ^^^d ^n.s{v) — 100) of temperature and polarization anisotropics 
(the normalization is chosen such that hn{rjr) = 1). 

one and the same mode function for both polarization states s, Eq. ()24() . Then, the coefficients 
a^(n, s), Eq. H76c|) . differ essentially only in sign for two different s, i.e. o£(n,-L) = —af_2{n,R). 
Therefore, their contributions will cancel out in expressions (jHEI) for Cf^ and Cf^. This statement 
is in agreement with Ref. \^\ . 

Without having access to a 'theory of everything' which could predict one unique distribution 
of the CMB radiation field over the sky, we have to rely on the calculated statistical averages (j86|l . 
We can also hope that our universe is a 'typical' one, so that the observed values of the correlation 
functions should not deviate too much from the statistical mean values. 

VI. EFFECTS OF RECOMBINATION ERA 

All our final graphs and physical conclusions in this paper are based on exact formulas and 
numerical calculations, starting from numerical representation of the key functions H(r]) and 
Eqs. (|55|) . (|6H) . However, we derive and explain all our results by developing manageable and 
accurate analytical approximations. At every level of calculations we compare exact numerical 
results with analytical ones. 

A. Temperature anisotropy angular power spectrum 

The temperature anisotropy power spectrum Cj^ is determined by the multipole coefficients 
aj{n,s), Eqs. (|86|) . (|76aj) . (|77aj) . The typical graphs for the functions Hn^siv) ^^'^ ^n,s{'n) 
shown in Fig. [l] 

Since the visibility function g{r]) is a narrow function, a convenient analytical approximation 
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is the hmit of an instanteneous recombination. The function e~'^ is replaced by a step function 
changing from to 1 at = Tj^ec, e""^ = h{r] — Tj^ec), and the function g{7]) is replaced by a delta- 
function, g{r]) = 6{r] — "qdec)- In this limit, the contribution to aj{n,s) from the scattering term 
^n^siv) is proportional to Arjdec- It can be neglected in comparison with the contribution from the 
gravitational term Hn^siv)- The ratio of these contributions is of the order of nAr]rec, and it tends 
to zero in the limit of instanteneous recombination, Aij^ec 0. 
Neglecting the scattering term, we write 

Vr s 

aj{n,s) = 7V4vr(2^ + l) / dv ^m). (87) 

Vdec 

This integral can be taken by parts, 

VR ^ 
hn (.'ndec)Te{Cdec) + J dj] n hn (^)~^^ 

Vdec 

The remaining integral contains oscillating functions and its value is smaller, for sufficiently large 
n's, than the value of the integrated term. This is illustrated in Fig. O Therefore, we have 



aj(n,s) =7V47r(2£ + l) 



aj{n, S) = -7^471(2^ + 1) hn (%ec)7>(Cdec). (88) 

Finally, we put Eq. (|S5|) into Eq. and take into account the definition of the metric power 
spectrum (|13p. Then, we get 

Cr = W I ^h\n, Vdec)Ti{Cdec). (89) 

The projection factor Tf{Cdec) is given by Eq. ()75a() . Since the spherical Bessel functions 
reach maximum when the argument and the index are approximately equal, Cdec ~ ^, a partic- 
ular wavenumber n is predominantly projected onto the multipole d k. n: 

^ ~ Cdec = ^{'qR - Vdec) ~ n. (90) 

This can also be seen in Fig. [51 Thus, the oscillatory features of the metric power spectrum 
h'^{n,r]dec) in the n-space are fully responsible for the oscillatory features of the angular power 
spectrum £{£ + \)Cj'^ in the ^-space (We use this opportunity to correct a misprint in Fig. 2 
of Ref. 1]| : the plotted lines are functions C;, not l{l + l)Ci. For the early graphs of Ci see Ref. 

The g.w. metric power spectrum h'^{n,rjdec) for the case /3 = —2, and the function £{£ + 1)C'['^ 
caused by this spectrum, are shown in Fig. [Tfa) and Fig. [7|^b). The normalization of the metric 
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FIG. 5: The contributions to the power spectrum t{l + l)Cj"^ from an individual mode n. The solid line 
shows the exact result calculated according to (I87II . while the dashed line shows the approximation H88|l . 
The normalization has been chosen such that hn(r\r) = 1. 

power spectrum is such that the function 1(1 + l)Cj'^ at ^ = 2 is equal to 1326 /iK^ The 
interval ^ < 90 is generated by waves with n < 90. These waves did not enter the Hubble radius by 
the time r/^ec- Their amplitudes are approximately equal for all n's in this interval (compare with 
Fig. The gradual decrease of the angular power spectrum at larger ^'s is the reflection of the 
gradual decrease of power in shorter gravitational waves whose amplitudes have been adiabatically 
decreasing since the earlier times when the waves entered the Hubble radius. 

B. Polarization anisotropy angular power spectrum 

The decisive function for polarization calculations is ^n,s(j])- We have approximated this func- 
tion by ^n'ls, Eq. (|73() . and compared it with exact result in Fig. |5J For qualitative derivations it 
is useful to make further simplifications. 

Since g(r\) is a narrow function, the integral in Eq. (|73|) is effective only within a narrow interval 

s 

Ar/rfec- Assuming that the function d hn {v)/dri does not vary significantly within this interval, we 
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can take this function from under the integral 



(91) 



Clearly, the assumption that the function d hn (^)/'^?? is almost constant within the window 
^Vdec gets violated for sufficiently short waves. With 1)91^ . we expect degradation of accuracy for 
wavenumbers n approaching n*. This is illustrated by a dotted line in Fig. |5J Nevertheless, the 
approximation (|91() is robust for n < n*, and it reveals the importance of first derivatives of metric 
perturbations for evaluation of the CMB polarization. 

We now introduce the symbol P to denote either E ox B components of polarization. In terms 
of multipoles (n, s), the angular power spectrum is given by Eq. l|86|) . Putting (|9T|l into Eqs. 
(|77b|) . H77cl) and denoting by ^^(C) the respective projection functions, we get 



af in, s) 



7 V47r(2£ + 1) j c 



dr] 



I dhn jr]) 
10 dr] 



PiiO 







(92) 



Again referring to the peaked character of g{rf) and assuming that the combination {d hn 
{f]) / dr})Pi[C,) does not change significantly within the window A-qdea we take this combination 
from under the integral, 

/ s 

I dhn (ll) 



af(n,s) « 7V47r(2£ + l) D{ 



n] 



10 dr] 



-m) 



A. 



(93) 



T] — Tjre.c 



The two new factors in this expression, D{ri) and A, require clarification. 

The factor D[ri) compensates for gradual worsening of our approximation when the wavenumber 
n approaches n^,. For large n's, the functions under the integrals change sign within the window 
^^dec5 instead of being constant there. This leads to the decrease of the true value of the integral 
in comparison with the approximated one. The evaluation of this worsening suggests that it can 
be described by the damping factor 



D{n) 



1 + 



nAr]de 



We inserted this factor 'by hand' in Eq. ()93p. Additional arguments on this damping are given in 
Ref. 3 (see also H)- 

The factor A is the result of the remaining integration over r] in Eq. (El |3> 




dr]'e~T5-M^ \=}^ I dr] 



Vr 



36 



Since e rapidly changes from to 1 around recombination, the integrand (e" 



-(3/10)r(,7) 



-T{ri) 



) 



is nonzero only there, and A is expected to be of the order of Arjdec- To give a concrete example, 
we approximate g{r]) by a gaussian function 



div) = , — 

V2^(A7?rfec/2) 

Then, the quantity A can be found exactly, 

+00 



A 



7 



At?, 



dec 



dx 



— \ — erf , ^ , 
2 2 



exp 



jr] - r]decf 
'2(A%ec/2)2 



— \ — erf 
2 2 



X 

71 



0.96A7]dec- 



Clearly, factor A in Eq. H93|) demonstrates the fact that the CMB polarization is generated only 
during a short interval of time around recombination. 

Finally, substituting (|93|) into l|86jl and recalling l|14j) . we obtain the polarization angular power 
spectrum: 

2 



•' s=l,2 



s 

d h {n,r}) 




drj 


n — ^rec 



PiiCr 



(94) 



Similarly to the case of temperature anisotropies, the projection factors P^(Cdec) predominantly 
translate n into i according to Eq. ()9U() . The oscillatory features of the power spectrum of the first 
time-derivative of metric perturbations get translated into the oscillatory features of the power 
spectra for E and B components of polarization. This is illustrated in Fig.[7|^c) and Fig. [7|^d). The 
waves with nr?^ec did not enter the Hubble radius by ?? = r/^ec- They have no power in the 
spectrum of dh{n,ri)/dr] at n <C 90, and therefore there is no power in polarization at ^ ^ 90. On 
the other hand, the first gravitational peak at n ~ 90 gets reflected in the first polarization peak 
at£^ 90. 

There is certain difference, however, between the projection functions Ei[C,dec) and Bi{C,ciec)- 
This is shown in Fig. El The B-mode projections are more 'smeared' and their maxima are shifted 
to somewhat lower ^'s. This explains the visible difference between Cf^ and Cf^ in Fig. [7| (see also 
The polarization angular power spectra, plotted in Fig. [3 were found from exact numerical 
calculations. 



C. Temperature-Polarization cross correlation 

A very special role belongs to the TE cross-correlation spectrum. We will show below that 
the TE correlation at lower €s must be negative for gravitational waves and positive for density 
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FIG. 6: Relative contributions of an individual Fourier mode n — 100 to various multipoles i in the power 
spectra i{l + l)Cf for temperature and polarization. The normalization has been chosen arbitrarily but 
same in both the graphs. 




FIG. 7: The left panel shows (a) the power spectrum of temperature anisotropics (.{t + l)Cj"^ (in /iK^) 
generated by (b) the power spectrum of g.w. metric perturbations Ijl^l) . (3 = —2. The right panel shows (c) 
the power spectra of polarization anisotropics + \)Cf^ (solid line) and l{t + \)Cf^ (dashed line), panel 
(d) shows the power spectrum of the first time derivative of the same g.w. field, Eq. H15|l . 
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perturbations. This distinctive signature can turn out to be more valuable for identification of 
relic gravitational waves than the presence of the B polarization in the case of gravitational waves 
and its absence in the case of density perturbations. The expected TE signal from gravitational 
waves is about two orders of magnitude stronger than the BB signal, and it is much easier to 
measure. At lower ^'s, the contributions to TE from gravitational waves and density perturbations 
are comparable in absolute value, so the g.w. contribution is not a small effect. The total TE 
cross-correlation has already been measured at some level P, |3| • 

To find the TE power spectrum we have to use the product of aj^{n,s) and af^{n,s) in Eq. 
(jHHJ. For qualitative analysis we will operate with the approximate expressions (|88j) and (|9.S|) . 
Then, the TE correlation reads 

(s s 
A / ^dh{n,r]) s d h* {n,r]) 

(Te{(dec)EiiCdec)) . (95) 

For a given £, the projection factor {T£{Qec)E£{Cdec) ) peaks at n ^ i and is positive there. There- 
fore, the sign of Cf^ is determined by the sign of the term: 



V=Vde 



^ h* t \dh{n,vi) s dh*{n,r]) 



dh?'{n, if) 
dij 



(96) 

V=Vdec: 



V=Vdec 

The adiabatic decrease of the g.w. amplitude upon entering the Hubble radius is preceded by 

s 

the monotonic decrease of the g.w. mode function hn (v) ^ ^ function of ry. This behaviour 
is illustrated in Fig. ^ It is clear from the graph that for n < 100 the quantity (|96|) is negative, 
because the first derivative of hn [f]) is negative. Therefore, for ^ ^ 90 the correlation Cj^ must be 
negative. For larger £'s the TE correlation goes through zero, changes sign and oscillates reflecting 
the oscillations of the function H96|) in the n-space. An exact numerical graph for TE correlation 
caused by gravitational waves is depicted in Fig. |HJ The graph clearly shows how the sign and 
features of the spectrum ()96() get translated into the sign and features of the TE correlation. 

It is shown in Appendix El that in the case of density perturbations the TE correlation must 
be positive at lower €s. This is because the relevant metric perturbations associated with density 
perturbations are growing in time and therefore the first time-derivative of metric perturbations is 

aitive. Because of other contributions, the TE correlation is expected to change sign at £ ~ 70 
, The region of intermediate multipoles 15 < £ < 90 should be of a particular interest. On 
one hand, the multipoles £ > 15 are not affected by the reionization era and its uncertain details. 
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FIG. 8: The bottom panel shows the spectrum H96() of gravitational waves, whereas the top panel shows 
the angular power spectrum £{£ + \)Cj^ caused by these waves. The negative values of these functions are 
depicted by broken lines. 

except the overall suppression by e~^'^''"=»°". On the other hand, at ^ < 90 the g.w. contribution to 
TE is not much smaller numerically than the contribution from density perturbations. The lower 
multipoles ^ ^ 15 are affected by reionization, and we shall study reionization in the next Section. 



VII. EFFECTS OF REIONIZATION ERA 



The reionization of the intergalactic medium by first sources of light has occurred relatively 
late, at z ~ 30 — 7 (see, for example, 0,13]). In contrast to the recombination part of the visibility 
function g{ri), which is narrow and high, the reionization part of g{rj) is broad and much lower 
(see Fig. I15j) . For a crude qualitative analysis, one can still apply analytical formulas derived for 
recombination. One has to replace rjdec with ijreion and Ar/^ec with Aijreion- The waves that start 
entering the Hubble radius at rjreion, i-e. waves with wavenumbers nijreion ~ provide the most 
of power to the spectrum of first time derivative of metric perturbations. Therefore, these waves, 
with n ~ 12, produce a 'bump' in polarization spectra at the projected £ ~ 6. However, these 
wavenumbers are comparable with the wavenumber n**, n^:^:Ar]reion ~ 27r, that characterizes the 
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FIG. 9: The source function ^n{v) rcionization era for different wavenumbers. 

width of the visibihty function at reionization. The assumption n <^ n** is not well satisfied, the 
analytical approximation becomes crude, and one has to rely mostly on numerics for more accurate 
answers. 

The numerical solution to the integral equation (jHJ) in the reionization era is shown in Fig. IHl 
One can see that the damping effect is expected to commence from n > 20, as the polarization 
source function ^niv) begins to show an oscillatory behaviour. The projection relationships are 
also far away from the almost one-to-one correspondence (. k, n that was typical for recombination 
era. In Fig. we show the contributions of a given n to various €s in the polarization power 
spectra + 1)C^. One can see that a considerable portion of power from a given n is distributed 
over many lower-order .^'s. 

The total effect of the reionization era is shown in Fig. 1111 This numerical result was based on 
our simplified model of homogeneous reionization, as described in Appendix^ The reionization 
'bumps' at lower multipoles for E and B components of polarization are similar in shape and 
numerical value. 



41 



Ijl !-. 



.xIOZ 



0.5 



10 20 30 40 

multipole I 





n = 20 







10 20 30 40 

multipole I 



10 20 30 

multipole I 




10 20 30 

multipole I 



FIG. 10: The contributions to the angular power spectra from an individual mode n. The dashed line is for 
e{£ + l)Cf^ and the solid line is for i{i + l)Cf ^. 
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FIG. 11: The rcionization 'bump' in the polarization power spectra. The dashed lines show the spectra 
when the rcionization era is ignored. The solid lines include both eras. The darker lines are for _B-mode, 
and the lighter lines for i?-mode. 
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VIII. COMPARISON WITH AVAILABLE OBSERVATIONS; SIGNATURES OF RELIC 
GRAVITATIONAL WAVES 



The theory that we are using here is appHcable to any primordial spectral index n. The initial 
conditions for gravitational waves, Eqs. H29|) . H30|) . as well as analogous initial conditions for density 
perturbations, hold /3, and hence n, as a free constant parameter. The spectral index n can be 
larger, equal, or smaller than 1. However, from the theoretical point of view, the 'red' primordial 
spectra n < 1 (/? < —2) seem to be unacceptable, or at least questionable. If /? < —2, the mean 
square fluctuations of gravitational field, Eq. H12|). are power-law divergent in the limit of very 
small wavenumbers n. One could argue that the extrapolation of the primordial spectrum to the 
'infrared' region of very small n's is uncertain, and for some reason the shape of the spectrum 
bends in the 'infrared' region making the integral 1)12(1 convergent at the lower limit. We prefer 
not to hide behind this possibility. If the shape of the primordial spectrum is allowed to be 
varied, then practically anything in the CMB data can be explained by the properly adjusted 
primordial spectrum. Therefore, our theoretical preference (unless the data will enforce us to 
change this preference) is a constant primordial spectral index n > 1 (/3 > —2). Obviously, such 
primordial spectra entail no difficulty in the 'ultraviolet' region of very large n's, because such 
short-wavelength fluctuations (today's wavelength ~ 3 cm.) did not satisfy the requirements of 
superadiabatic amplification and simply have not been generated. 

We now return to the TE correlation. In Fig.Elwe show the contributions to the TE correlation 
function from relic gravitational waves (gw) and primordial density perturbations (dp). (In order 
to enhance lowest ^'s we use the combination {i + l)C£ rather than i{£ + l)Ci.) For this illustration, 
we choose a flat primordial spectrum n = 1 (/? = — 2) and assume equal contributions from (gw) 
and (dp) to the temperature quadrupole: R = 1, where 

R^CfJ,i9w)/CjJ,{dp). 

We include the effects of reionization according to the model with Treion = 0.09. The (gw) contri- 
bution is numerically calculated from the solution to the integral equation ()61() . whereas the (dp) 
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3l|. One can see from the graph that 



contribution is plotted according to the CMBfast code 
the negative TE correlation function at lower ^'s is only possible if there is a significant amount of 
primordial gravitational waves. One can also see from the graph that a mis-interpretation of the to- 
tal TE eff^ect as being caused by density perturbations alone, could lead to a serious mis-estimation 
of 
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451 shows clusters of 



The TE correlation at lower £'s measured by the WMAP mission 
data points, including the negative ones, that lie systematically below the theoretical curve based 
on density perturbations alone. It is true that the data points in the interval 10 < ^ < 70 are 
concentrated near a zero level, the error bars are still large, and the measured TE correlation can 
be appreciably different from the theoretical statistically averaged TE correlation. However, the 
recent paper page 26, explicitly emphasizes the detection of the TE anticorrelation by WMAP: 
"The detection of the TE anticorrelation near ^ ~ 30 is a fundamental measurement of the physics 
of the formation of cosmological perturbations...". As we have already stated several times, the 
TE anticorrelation at lower ts, such as ^ ~ 30, can only take place (within the framework of all 
other common assumptions) when a significant amount of relic gravitational waves is present. 

Our theoretical position, as explained in this and previous papers, is such that we are asking 
not 'if relic gravitational waves exist, but 'where' they are hiding in the presently available data. 
We shall now discuss some models that fit the CMB data and contain significant amounts of relic 
gravitational waves. More accurate observations with WMAP, and especially Planck, should firmly 
settle on the issue of the sign and value of the TE correlation at lower multipoles. Hopefully, these 
observations will establish the presence of relic gravitational waves beyond reasonable doubts. 

To sharpen the discussion of allowed parameters, WG take the model with, n — 1.2, '^reion — 

0.09 

and R = l. We take the values ^mh? = 0.03, l^^/i^ = 0.12, h = 0.75, and we normalize the g.w. 
contribution to ^^tii = 440 uK^ at ^ = 2. In Fig. [T3k we show our calculation of the TT 

Q n 

correlation function in comparison with WMAP data and the best fit ACDM model [1^. (If it 
comes to the necessity of explaining 'dark energy', natural modifications of general relativity will 
be superior to unnatural modifications of the matter sector [s^.) One can see from the graph that 
even this model (which lies, arguably, on a somewhat extreme end) is consistent with the TT data 
at all ^'s and significantly alleviates the much discussed tension between theory and experiment at 
^ = 2. In Fig. 113b we show the TE correlation for exactly the same model. One can see that the 
inclusion of relic gravitational waves makes more plausible the negative data points at lower €s. 
Since the relative contribution of gravitational waves becomes small at ^ > 90, the higher £ portion 
of the graph is governed by density perturbations alone. Obviously, models with a little smaller n 
or R reach the same goals. 

Finally, in Fig. we combine together all correlation functions induced by relic gravitational 
waves. The graphs are based on the discussed model with n = 1.2, i? = 1. The future detection 
of the BB correlation will probably be the cleanest proof of presence of relic gravitational waves. 
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FIG. 12: The dotted line shows the contribution of density perturbations alone, and the dashed line shows 
the contribution of gravitational waves alone. The solid line is the sum of these contributions. It is seen 
from the graph that the inclusion of g.w. makes the total curve to be always below the d.p. curve. 
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FIG. 13: The dashed line shows the best fit ACDM model without gravitational waves. The solid line shows 
a model with spectral index n = 1.2 and gravitational waves R = 1. 



0.1 /iK^ in the 



For the discussed model, the predicted level of the BB signal is ^^'^p-Cf^ 

region near the first peak at ^ ~ 90. This level of i3-mode polarization shouiO De aetectaDie Dy 
the experiments currently in preparation, such as CLOVER 51], BICEP 52] and others (see also 
[^). Obviously, in this paper, we ignore many complications, including astrophysical foregrounds. 
More information can be found in Ref. 



45 




FIG. 14: The summary of CMB temperature and polarization anisotropics due to relic gravitational waves 
with n = 1.2 and R = 1. 

IX. CONCLUSIONS 

In this paper, we summarized the properties of relic gravitational waves that are directly rele- 
vant to CMB calculations. We explained the reasons why we are working with g.w. backgrounds 
possessing specific statistical properties, amplitudes, spectral indices, etc. Then, we worked out, 
essentially from first principles, a theory of CMB anisotropies induced by relic gravitational waves. 
Some parts of this theory are rederivations and confirmations of previous studies, some parts are 
new. The important advantage of our approach, as we see it, is a transparent physical picture. We 
believe we have demonstrated in the paper that every detail of the derived correlation functions is 
under full analytical control. Clear understanding of the participating physical processes has led 
us to the conclusion that the TE correlation in CMB can be a valuable probe of relic gravitational 
waves. We compared our theoretical findings with the WMAP data. We believe that the TE 
anticorrelation detected by WMAP at ^ ~ 30 is certain evidence for relic gravitational waves in 
the already available data. We propose more accurate observations of the TE correlation at lower 
ts and believe that these observations have the potential of providing a firm positive result. 
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APPENDIX A: POLARIZATION STATES AND RANDOMNESS OF GRAVITATIONAL 
WAVES 



As stated above, the quantities | hn (vn describe the magnitude of the mean-square fluctuations 
of the g.w. field in the corresponding polarization states s. We consider a particular mode n of the 
field: 

/iij(n,r/,x) = Pij (n) U„ (??)e*"-'' +hn (f?)e~*"''' c„ j . (Al) 

s=l,2 ^ ^ 

The lowest-order independent g.w. correlation functions amount to 

{hij{n,rj,^)p'^{n) /i^^n', 7?, x') P '^'(n')) = 8 |L {r])f 5,s'S^'\n - n'). (A2) 

The action of the random g.w. field (jAlj) on free particles leads to their relative oscillatory 
motion. We refer this motion to a local inertial frame. Let a ring of free particles to lie in the 
(1, m)-plane; the ring encircles the axis n. Then, the mean-square amplitude of oscillations in the 
'+' polarisation state is determined by | hn P, whereas the mean-square amplitude of oscillations in 

2 

the 'x' polarisation state is determined by | /in p. In general, the random gravitational- wave field 

12 12 

can be such that the oscillation amplitudes | /in P and | /in P are different. But if | /in P / I /in 
then the (averaged) observed picture of oscillations is not symmetric with respect to rotations 
around n. 

Formally, the correlation functions of the field with | /in | 7^ | /in | do not have symmetry with 
respect to the change of polarisation basis. Indeed, the transition to the primed basis according to 
Eq. brings the gravitational wave mode (jAll) to the form 

/i,,(n,r/,x) = Y P'^J (n) (bn (^?)e*"-" + L*(r/)e-*"-^) , 

s=l,2 ^ ^ 

where 

1 1 1 2 2 

bn iv) = cos2ii hn {ri) Cn + sin 2V' hn iv) Cn, 

2 1 1 2 2 

bn iv) = - sin hn (i?) Cn + COS 2ip hn iv) Cn • 
Taking into account the relationships (jlOjl . we can now derive the correlation functions for the 
new polarisation components: 

{bnb*n') = (^\ hn l"^ cos^ 2tP + \ hn P sin^ 2V'^ 6^^^n-n'), 



{bj*n') = (^1 hn psin2 2V' + | L pcos2 2v) S'-^^n-n'), 

1 2 / 1 2 \ 

{bnb*n') = - I /in P - I /in M sin 2-0 COS 2^^ 5(3)(n-n'). 



(A3) 



47 



It is seen from (IA3|) that, in general, the ■0-dependence survives, and the assumption of statistical 
independence of polarisation components in one basis is not equivalent to this assumption in another 
basis. However, one recovers the original correlation functions <\A2^ from (|A3|) if the conditions 
((TH|) are fulfilled. 

Similar properties hold true for circular polarizations. A g.w. mode /ijj(n, 7y,x) expanded over 
circular polarisation states is given by 

/iij(n,ry,x) = hij{n,r],x) + /i*(n, 77, x), 



where 



hij{n,7],x)= (n) L (??)e*"''' Cn . 



(A4) 



=L,R 



We assume that the complex random coefficients Cy^ (s = L,R), satisfy the statistical conditions 
(dEI). 

The relevant independent correlation functions are calculated to be 

(/ii,(n,?7,x)p*^*(n) /i^,(n',7?,x')p''(n'))=4 L (v)^ ^s.'^^'Hn - n'), {s = L,R). (A5) 



If the observer views the motion of test particles from the — n direction, i.e. against the direction 

R ^ 
hn (r/) 



of the incoming gravitational wave, the function 



is responsible for the mean-square am- 



L 



hn iv) 



IS 



plitude of the right-handed (clockwise) rotations of individual particles. The function 

responsible for the left-handed (anti-clockwise) rotations. (For more details about the motion of 

free particles in the field of gravitational waves, see ^^■) 

Expansion HA4() preserves its form under transformations if one makes the replacements: 

h 'n =hn e^'^^ , h 'n =hn c"*^'^. Therefore, the correlation functions HA5|) do not change, regardless 

2 



the value of the amplitudes 



hn iv) — L,R). On the other hand, discrete transformations H22|l 



L R R L 

generate the replacements: h n =hn, h n =hn- Therefore, the sense of correlation functions (|A5|1 
changes from L to R and viceversa. The symmetry between left and right is violated, unless the 
the conditions (1231) are fulfilled. 



APPENDIX B: ASTROPHYSICAL PREREQUISITES 

(i) Ionization History 

The ionization history of the Universe enters our equations through the density of free electrons 
Ne{r]), i.e. electrons available for Thompson scattering. Specifically, we operate with the quantity 
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q{7]): 

g(r?) = aTa{r])Neirj)- (Bl) 
The optical depth r between some instant of time rj' and a later instant rj is defined by the integral 

V 

r(r?,r?')= / dr,"q{v"). 



rj' 



The optical depth from some rj to the present time r]ji is denoted r(r/) and is given by 

Vr 

Tiv) = t{vr,v) = / dr]'q{7]'). 



V 

It follows from the above definitions that r(r/, rj') = T{rj') — T(rj). 

A key role in our discussion is played by the quantity g{rj) called the visibility function: 

9{V) = q{ri)e-^^^'^ = -e'^^'^ , / g{ri)drj = 1 
dv Jo 

The state of ionization is determined by microphysical processes during all the evolution of 



the Universe 



5£, 



59|. For our purposes it is sufficient to focus on two eras: early recombination 
and late reionization. The recombination of primordial plasma into atomic hydrogen and helium 
is accompanied by decoupling of CMB radiation from the rest of matter. (For simplicity, we do 
not distinguish here the notions of recombination and decoupling.) This relatively quick process 
has happened at redshifts around z^^c ~ 1100. Much later, at redshifts around Zreion ~ 

10 9 the 

intergalactic medium has become ionized again, presumably, by ionizing radiation of first condensed 
objects. 



The density of free electrons is modelled 

Ne{ri) = (l 



60| by the expression 



2 J nip \ a{ri) 

where Yp ~ 0.23 is the primordial helium mass fraction, Xe{i]) is the fraction of ionized electrons, 
0,b is the baryon content, and nip is the mass of a proton. In the framework of linear perturbation 
theory it is sufficient to regard the electron density as homogeneous, i.e. depending only on ry. For 
Xe(?y) we use the fitting formula jo^: 

Yr,\~^ / Co \ / nir, \ „^., _1 / Z \C2-1 / a' 



^'^ V 2y \WOOJ \2aTlHPcJ ^ Viooo/ V« ^ ^ ' ^ ' 

where ci = 0.43, C2 = 16 + 1.8lnO,B, and z is a redshift. 
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FIG. 15: The graphs (a), (b), (c) show the quantities q{ri), T{ri), g{r}). The graph (d) is a zooming oi g{rj) 
in the region of recombination. The dashed hne shows a model without reionization. The solid line takes 
rcionization into account, with Treion = 0.1. 

As for the reionization, we assume that it was practically instantaneous and happened at a 
redshift Zj-eion ~ 16. The function X(. is determined by Eq. (|B2|) for z > Zreion and = 1 for 

Z ^ Zreiorf 

To plot the graphs for q{ri), t{i]) and g{ri) in Fig. 1151 we use the above-mentioned parameters 
and rjf, = 0.046 P, It is seen from Fig. that the visibility function g{7]) is sharply peaked 
at the era of recombination. The peak can be characterized by the position of its maximum r]dec, 
and the characteristic width Arjdec- A similar, but less pronounced, peak is present also at the era 
of reionization. 

(ii) Frequency dependence of the Stokes parameters 

As is seen from Eq. ()46() . the frequency dependence of both, temperature and polarization, 
anisotropies is governed by the function f{i)) 3i]. We assume that the unperturbed radiation field 
has a black-body spectrum. 

Ml") = 7 — ^— , 

exp I hu I kT I — 1 
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where T = Ta{rj) and the present-day value of T is T{r]ff) ~ 2.73K. 

It follows from these formulae that /(P) is approximately 1 in the Reyleigh- Jeans part of the 
spectrum, and /{u) varies as hu/kT in the Wein part. In practice, we are usually interested in 
the total intensity, i.e. Stokes parameters integrated over all photon frequencies. This integration 
produces the factor 7: 



7 



dv -^nQ{u)f{u) = -4/0, 



(B3) 



where /q = / ^^^n^iy) is the total intensity of the unperturbed radiation field. The factor 7 
often appears in the main text of the paper. 

APPENDIX C: TWO ESSENTIAL VARIABLES FOR TEMPERATURE AND POLAR- 
IZATION 



It is seen from Eq. H48|) that the frequency dependence of n'-"'"^ is determined by the factor 
/(t')no(t'). Therefore, we can single out this factor and write 

n^^^(??,z^,/i,(/>) = ^/(z>)no(i>)n(^)(??,/x,(/)). 

The equation that follows from 1)48(1 and (|49() reads 
d 



dr] Att 



+ 1 2tt 

^"^^ f f'dll'dcP'PilI,cp-f,',cP')h(^\7^,f,',<P'). 



(CI) 



-1 



The quantity n^-*-) (rj, (p) as a function of 4> can be expanded in a series 



+00 



m=—oo 



The explicit structure of the Chandrasekhar matrix P(^, /u', is given by the expression 
(equation(220) on p. 42): 

P(^, /i', 0') = Q [PW + (1 - ^2)^ (1 _ ^'2)ip(l) + p(2) 

where matrices Q, p(°)(/i; /i'), P^^\f^,4>; fi', (/)'), P(^)(^, /i', 0') read 

/l o\ 



(C2) 

Q 



Q 



1 
y 2y 
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p{0) 



/2 



1 




^ Afifi' cos I 



p(i) 



sm I 



— 2/ism (</)' — (f 



COS {(j)' - (/)) 



J 



p(2) = 2. 
4 



^ /iV^ COS 2 ((/.'- 0) (,Qs 2 - 0) -^2^' sin 2 ((/.'- ^ 

-^'2 cos 2 (0' - 0) cos 2 (0' - 0) ^' sin 2 {(j)' - 4>) 
y /i;u'2 sin 2 (0' — (/)) —fi sin 2 — 0) ^/x' cos 2 — (p) J 



The structure of the Chandrasekhar matrix P(^, (/>; /i', 0') is such that it does not mix the 
m(/)-dependence, that is, 

+ 1 2tt 



1 



I\„im4> ^imd 
~ e 



-1 



Moreover, this integral vanishes for all \m\ > 2. This means that Eq. ((HH) is a homogeneous 
differential equation for all with m ^ ±2. Assuming zero initial conditions at some initial rj, 
we obtain fi^^ = for all m / ±2. Hence, we are left with three functions of {r],iJ,): 

nW(r?,/i,<A)=nW(r?,/i)e±2^f (C3) 

We are now able to show that only two out of the three functions h^^\r],fj,) are independent. 
Indeed, using (|C3|) in Eq. (|C1() we arrive at a system of three linear equations for the components 
of n(-'-)(r/, /i): 



^ + qiv) + 



ni(?7,/x) 



\ 



1 



3 , ^ 



( ^2 \ 



J(7?), (C4) 



Y ±2i/i y 

where T{r\) is the remaining integral that was left over from the last term in Eq. HClf) : 



'^W) = o / t^At' ^lij], - h2{ri, fi') ± ifi'nsir], fi') 



(C5) 
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Despite the complicated appearance, only two of the three equations HC4|) are really coupled. 
Indeed, making linear combinations of equations HC4() it is easy to show that the combination 
2ifi{hi — n2) =F (1 + A'^)^3 satisfies a homogeneous differential equation. Assuming zero initial 
conditions, we derive 

2ifi{hi - 712) =F (1 + = 0. 



As two independent and essential variables we choose 



In terms of a and /? equations of radiative transfer reduce to Eq. (|51j) and Eq. I|52|). and the 
definition (|(]5|) for I{r]) takes the form of Eq. (|53|) . 

Although we have considered a gravitational wave perturbation, the existence of only two es- 
sential variables is a general statement and it applies to density and rotational perturbations as 
well. In general, the perturbed radiative transfer equation contains an arbitrary function /(r/, /i, (p) 
in front of u, rather than a specific combination (1 — ]i'^)e^'^^'^dhn/ drj quoted in Eq. HC1|1 . Function 

/(ry, /X, 4>) can be expanded in a series similar to Eq. (|C2|1 : 

+00 

m=—oo 

Since the scattering integral (second term in the r.h.s of Eq. 1)011) ) vanishes for all \m\ > 2, the 
functions hm\'r],fJ') are fully determined by fmiv^t^) > 2) and describe the temperature 

variations only. To discuss polarization, we have to consider three remaining cases m = 0, ±1, ±2. 
The explicit identification of the two essential variables a{r], fi) and (3{rj, n) for the case m = ±2 

L 

has been given above. Specifically for gravitational waves, /2 = (1 — l^ )dhn f-2 = 

R 

(1 - fi^)d hn {r])/dr]. 

The identification of the two essential variables for the case m = proceeds in a similar manner. 
Having calculated the scattering integral for m = one can show that the equation, analogous to 
Eq. (|C4() . will now read 



d_ 
where 



\ 



n2(??,//) 



h{Ti,ij) + <i{v)io\'n) 



1 



3 / X 



3;u2 - 2 
1 




\ 



+1 
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+1 



Jiv) = ^ I dfJ-' - 2)ni(r/,/i') +ri2(r?,/x') 



-1 

Obviously, ns satisfies a homogeneous equation and can be put to zero. Specifically for density 
perturbations (see Appendix function /o(r/,/i) consists of terms representing gravitational field 
perturbations and the Doppler term arising due to the baryon velocity: 

We have checked that in the case m = ±1 the problem also reduces to only two essential 
variables. The combination ni — n2 =F ifJ-ns vanishes at the zero initial data. 

APPENDIX D: TEMPERATURE AND POLARIZATION ANISOTROPIES CAUSED 
BY DENSITY PERTURBATIONS 

To discover relic gravitational waves in the CMB data we have to distinguish their effects 
from the effects of density perturbations. The theory of temperature and polarization anisotropies 
caused by primordial density perturbations is very much similar to the theory of relic gravitational 
waves. We start from the metric Fourier expansion (|10)) with the polarisation tensors Q. Having 
derived and solved integral equations of radiative transfer in the presence of density perturbations, 
we arrive at our final goal of distinguishing the TE cross-correlations. 

(i) Radiative transfer equations 

The equations of radiative transfer in the presence of a single mode n of density perturbations 
are similar to Eq. ()48() and read 



9 / N . 
^ + qiv) + Ui 
or) 



-(1)/ - i 



/(z>)no(i>) 



(n) — -j^ g(r?)e Vi 



Qiv) 
47r 



(Dl) 



where the extra term e^Vi takes care of the movement of scattering electrons with respect to the 
chosen synchronous coordinate system 

For technical reasons, it is convenient to work with the 'scalar' h{ri) and 'longitudinal' hi{ri) 
polarization mode functions, instead of the original hn {fj) 3l- The relationship between them is 



hn iv) 



l{h{rj)-hi{r,)), hniv) 
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where the wavenumber index n on h{r]) and is imphcit. Both polarization components of 

metric perturbations participate in Eq. (|D1|) . In the frame associated with the density wave, i.e. 

■ ■ * __ 

for n/n = (0,0, 1), the structures e^e^ Pij (see Eq. lOJ) and e'^Vi = —ifJ^vi, depend on /x = cos 6, but 
not on the azimuthal angle <j). 

By argumentation similar to that in Appendix O one can show that a solution (not vanishing 
on zero initial data) to Eq. (jPlj) must have the form 

"n(?y,^) 



f{i')no{v) 



1 



V / 



(D2) 



Substituting Eq. ()D2|) into Eq. HDljl we arrive at a system of coupled equations for a and /5 



d 



drj 



a„(?7,/i) 
/?n(??,At) 



1 /d/i 



2 Vd?/ 
1 



1 



g(r/) 1- P2{^i)]l2 



(D3a) 
(D3b) 



where 



+1 



^liv) = ^ I dfi an{r],fj,), 



-1 
+1 



^2{v) = 2 I ^1-^ (1 -^2(/^)) /?«(?/, At) - P2{fJ-)an{ri,fi) 



(D4a) 



(D4b) 



The quantity Ti is the monopole component of the perturbed radiation field, whereas T2 is the 
quadrupole {£ = 2) component, responsible for the generation of polarization. 

To make contact with previous work we note that the variables a and f3 are closely related to 
the variables At and Ap from Ref. [s]]. Assuming a black body unperturbed radiation field, we 

d. We also note that the mode functions h and hi are related to the mode 

nn 

,|6J| hy h = —2r], /i; = — (h + 67]). Keeping in mind the difference in 
notations, one can verify that equations (|D3a|) and HD3b|) are equivalent to equations (11) in Ref. 

Q. 

(ii) Integral equations and their solutions 



have At = —a, Ap - 
functions h and r] used in 
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A formal solution to equations (|D3a|) and (|D3b|) can be written as 
an{'n,lj) = J dr]'e 





n' p-'riV,-n')--i-ntJ.{ri-ri') 



(D5a) 
(D5b) 



Proceeding in a manner similar to that in Sec. IIVI we substitute (|D5a|) and HD5b|) into HD4a|) and 
(|D4b|) . After certain rearrangements we arrive at two coupled integral equations for Zi and X2: 

V 







1 



^i{v')jo{x) + -l2{ri')j2{x 



V 



+ J dr] e 




1 f dh dhi d^jo > , M / '\ - f \ 



(D6a) 



^2{v) = J dr]'q{r]')e 




+ J drj e 




1 f dh , . dh. , , , 

2 W''"'^*7 1 



2x" 



« / ■ / N , ^(^^kW , I IX f ixdj2 

+ q{v )Vb{v )■ 



dx^ 



dx 



(D6b) 



where the argument of spherical Bessel functions je{x) is x = n{r] — r]'). 

Equations (|D6a() and HD6b|) together with the continuity equations for matter perturbations 



and Einstein equations for metric perturbations (see for example 
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6l 



63, 



65 [ 



, form a closed system 
, the radiative transfer 



of coupled integro-differential equations. In previous treatments 
equations were presented as an infinite series of coupled ordinary differential equations. 

Similarly to what was done in Sec. lIVBl we can analyze Eqs. ()D6a|) . HD6b|) in terms of expansions 
in powers of n, 



k=0 



T2 = J^xfV^ 

k=0 



The kernels of Eqs. ()D6a|) . l|D6b|) are expanded in powers of x. In the long- wavelength approx- 
imation, i.e. for x ^ 1, we can limit ourselves by first non-vanishing terms. The zero order 
approximation brings us to the equations 
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-— + 2q(r] jvbir] )x 



The solution to these equations is given by 



4'\v) = J dv' 





4"/*'"'<'''>-'"''"'-c"' 





(D8a) 



15 



dhi 
drj 



1 + '^^ j d7]"q{r]")e 



(D8b) 



(iii) Multipole coefficients 

We are mostly interested in the present-day values of a and (3 and, hence, we put ?? = r//j in 
HD5a|) and ()D5b|) . Irrespective of approximation in which the functions ^2{fl) ^^'^ known, 

the multipole coefficients for the radiation field can be found in a way similar to that in Sec. 



«L(^) = {-iYl^/^^{2l + l)6^oaJ{n), 
^fmi^) = (-i)'7\/4^(2^ + l)'5™oaf(n). 



af^{n) = 0, 



where aj{n), af{n) are given by 

Idh 1 dhi d^ 



VR 



aj{n) 



afin) 



drj 



2)! 



1 »?« 



(D9a) 
(D9b) 



and C = n(77B — r?). As expected, in the case of density perturbations all = for m 7^ and 
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aL(n) = 



l32l.l67j. The formal reason for this is that all quantities in HD2j) do not depend 
on (f) (A more detailed exposition can be found in [^, Q). 

Let us start from the temperature multipoles. For illustration, we consider an instanteneous 
recombination. We replace e""^ with the step function, g{rj) with the (5-function, and we neglect 



the I2 term. Then, we get 



aj{n) 



2^1 + ^^6^ ) hiO 



Vr 



+ 



V=Vdec 



drj 



fdh dhi d^ 
\dri dr] dC,"^ 



Vde 
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This expression can be further simphfied by taking the remaining integral by parts. After some 
rearrangements we arrive at the final expression 



aj{n) 



+ 



2v? drf' 6re drj 

1 d'^hi 



2n^ drj^ 



+ Xi + 



1 dhi 
2n dr] 



+ Vb 



dC 



MO 



V=11dec 



VR 



drj 



drj drj^ ' "^^ 



(DIO) 



ride 



For the growing mode of density perturbations in the matter-dominated era, the integrand of 
the remaining integral vanishes. If, in addition, the intrinsic temperature perturbation Ii and 
the plasma velocity Vh are zero at r? = rjdec-, we recover from Eq. (|Dlfl|) the four terms of the full 
Sachs- Wolfe formula (43) in Ref. _68] . We would like to note in passing that the origin of the often 
used, including this paper, combination + l)Ci is in fact a historical accident, when a wrong 
motivation leads to a convenient notation. This combination arises in the essentially incorrect 
formula + 1)C^ = const that can be derived after the unjustified neglect of three out of four 
terms in the original Sachs- Wolfe formula (43) in Ref. ^^I- This clarification is important not only 
by itself, but also for correct physical interpretation of reasons for the rise toward the 'first peak' 
at ^ ~ 220 of the actually observed function + 1)C^. For more details, see jod ]. 

We now turn to polarization anisotropies. We use the zero order approximation ()D8b|l in HD9b|) . 

/ 



1 ({l + 2)\ 



20 \{l-2)\ 



1 -riR 
2 



dr] ge^^"^ / dr]'e 



\ 



dhi 
drj' 



j + 2n dr]"q{r]")e 



-r{v',v") 



Vbiri") 



As g{r]) is a narrow function, the integrand is localized near r] = %ec- By the same argumentation 
that was used in Sec. lVI Bl the above expression can be reduced to quantities evaluated at r] = rjdec-, 



af{n) 



{£ + 2)1 
20 V(^-2)! 



MO 
C2 



dhi , r\ 
dr] J 



(Dll) 



V=Vdec 



where the factors A and A are of the order of the width of g{rj) in the recombination era. Explicitly, 
they are given by 



vr 



VR 



A 



dll giv) / dr]' e w"^ 



A 



dr] g{r]) j drj' e-To^('''^') j drj" g(7?")e-^(^''^"). 


The angular power spectra in the case of density perturbations are described by formulas anal- 
ogous to ()86() . By the steps similar to those in Section fV CI it can be shown that the power spectra 
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are given by 

(jxx' ^ ^ / ndnaf{n)af'{n), (D12) 



vr 

where C for density perturbations is C = \/247rZp/. 
(iiii) Temperature-polarization correlations 

We focus on the TE correlation at the relatively small multipoles, H. < 70. The dominant 
contribution to these multipoles comes from density perturbations that did not enter the Hubble 
radius at recombination, n < 70. It is sufficient to consider the early time evolution of these 
perturbations in the matter-dominated era. 

Restricting our analysis to the growing mode, we write for the metric and matter perturbations 
(see, for example, S^): 

ST 2 2 
h{r]) = Bn = const, hi{ri) = Bn{nrjf, z[^^ = — = - -Bn{nr]f , Vb = - -Bn{nr]f . 



Substituting this solution into (|D10|) . (|D1H1 and taking into account only the lowest-order terms 
in n (the long-wavelength approximation) we arrive at 

aj ~ -Bn MO , (D13a) 

V=Vdec 



Finally, from Eq. (|D12() we obtain for the TE correlation: 



(D13b) 



^ 10 V(^-2)!y 7 " V C 



(D14) 



V=Vdec 

All terms in the above expression are strictly positive. Therefore, the TE correlation caused 
by density perturbations, in contrast to the case of gravitational waves, must be positive at lower 
multipoles ^ ^ 50 where our approximations are still valid. As mentioned before, this difference 
between gravitational waves and density perturbations boils down to the difference in the sign of 
first time-derivative of the associated metric perturbations. 
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Abstract 

' A strong variable gravitational field of the very early Universe inevitably generates relic gravitational waves 

_ by amplifying their zero-point quantum oscillations. We begin our discussion by contrasting the concepts 

O I of relic gravitational waves and inflationary 'tensor modes'. We explain and summarize the properties 

^ ■ of relic gravitational waves that are needed to derive their effects on CMB temperature and polarization 

anisotropics. The radiation field is characterized by four invariants /, V, E, B. We reduce the radiative 
transfer equations to a single integral equation of Voltairre type and solve it analytically and numerically. 
, We formulate the correlation functions Cf^ for X, X' ~ T,E,B and derive their amplitudes, shapes and 



oscillatory features. Although all of our main conclusions are supported by exact numerical calculations, we 
obtain them, in effect, analytically by developing and using accurate approximations. We show that the TE 
correlation at lower i"s must be negative (i.e. an anticorrelation) , if it is caused by gravitational waves, and 
positive if it is caused by density perturbations. This difference in TE correlation may be a signature more 
valuable observationally than the lack or presence of the BE correlation, since the TE signal is about 100 
times stronger than the expected BB signal. We discuss the detection by WMAP of the TE anticorrelation 
at £ ~ 30 and show that such an anticorrelation is possible only in the presence of a significant amount 
of relic gravitational waves (within the framework of all other common assumptions). We propose models 
containing considerable amounts of relic gravitational waves that are consistent with the measured TT, TE 
and EE correlations. 

PACS numbers: 98.70.Vc, 98.80.Cq, 04.30.-w 
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I. INTRODUCTION 



The detection of primordial gravitational waves is rightly considered a highest priority task for 
the upcoming observational missions 1]. Relic gravitational waves are inevitably generated by 
strong variable gravitational field of the very early Universe. The generating mechanism is the 

n 

superadiabatic (parametric) amplification of the waves' zero-point quantum oscillations ;2]. In 
contrast to other known massless particles, the coupling of gravitational waves to the 'external' 
gravitational field is such that they could be amplified or generated in a homogeneous isotropic 
FLRW universe. This conclusion, at the time of its formulation, was on a collision course with the 
dominating theoretical doctrine. At that time, it was believed that the gravitational waves could 
not be generated in a FLRW universe, and the possibility of their generation required the early 
Universe to be strongly anisotropic (see, for example, P]). 

The generating mechanism itself relies only on the validity of general relativity and quantum 
mechanics. But the amount and spectral content of relic gravitational waves depend on a specific 
evolution of the cosmological scale factor (classical 'pumping' gravitational field) a{rj). The theory 
was applied to a variety of a(rj), including those that are now called infiationary models (for a 
sample of possible spectra of relic gravitational waves, see Fig. 4 in Ref. 4]). If a unique a{r]) 
were known in advance from some fundamental 'theory of everything', we would have derived the 
properties of the today's signal with no ambiguity. In the absence of such a theory, we have to 
use the available partial information in order to reduce the number of options. This allows us to 
evaluate the level of the expected signals in various frequency intervals. The prize is very high - 
the actual detection of a particular background of relic gravitational waves will provide us with 
the unique clue to the 'birth' of the Universe and its very early dynamical behaviour. 

A crucial assumption that we make in this and previous studies is that the observed large- 
scale cosmic microwave background (CMB) anisotropies are caused by cosmological perturbations 
of quantum-mechanical origin. If this is true, then general relativity and quantum mechanics 
tell us that relic gravitational waves should be a significant, if not a dominant, contributor to 
the observed large-scale anisotropies. From the existing data on the amplitude and spectrum of 
the CMB fiuctuations we infer the amplitude and spectral slope of the very long relic gravitational 
waves. We then derive detailed predictions for indirect and direct observations of relic gravitational 
waves in various frequency bands. 

At this point, it is important to clarify the difference between the concepts of relic gravitational 
waves and what is now called inflationary gravitational waves. The statements about infiationary 
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gravitational waves ('tensor modes') are based on the inflation theory. This theory assumes that 
the evolution of the very early Universe was driven by a scalar field, coupled to gravity in a special 
manner. The theory does not deny the correctness of the previously performed calculations for 
relic gravitational waves. However, the infiationary theory proposes its own way of calculating the 
generation of density perturbations ('scalar modes'). The infiationary theory appeals exactly to 
the same mechanism of superadiabatic (parametric) amplification of quantum vacuum fluctuations, 
that is responsible for the generation of relic gravitational waves, but enforces very peculiar initial 
conditions in the 'scalar modes' calculations. 

According to the inflationary initial conditions, the amplitudes of the 'gauge-invariant' metric 
perturbations (" associated with the 'scalar modes' (or, in other words, the amplitudes of the 
curvature perturbations called ( or TZ) can be arbitrarily large from the very beginning of their 
evolution. Moreover, the theory demands that these amplitudes must be infinitely large in the 
limit of the deSitter expansion law a(r/) oc which is responsible for the generation of a fiat 

(Harrison-Zeldovich-Peebles, 'scale-invariant') primordial spectrum, with the spectral index n = 1. 
At the same time, the amplitudes of the generated gravitational waves are finite and small for 
all spectral indices, including n = 1. Since both, gravitational waves and density perturbations, 
produce CMB anisotropies and we see them small today, the infiationary theory substitutes (for 
'consistency') its prediction of infinitely large amplitudes of density perturbations, in the limit 
n = 1, by the claim that it is the amount of primordial gravitational waves, expressed in terms 
of the 'tensor /scalar ratio r', that should be zero, r = 0. For a detailed critical analysis of the 
inflationary conclusions, see P]; for arguments aimed at defending those conclusions, see Q]. 

The science motivations and CMB data analysis pipelines, designed to evaluate the gravitational 
wave contribution, are usually based on inflationary formulas ^jQjyj^^l- In particular, according 
to the inflation theory, the primordial power spectrum of density perturbations has the form (it 
follows, for example, from Eqs. (2.12a) and (2.12b) in Ref.|lo| or from Eqs. (18) and (19) (or 
(A12) and (A13)) in Ref.Q): 

Despite the fact that this spectrum diverges at r = (n^ = 1, nj = 0) the CMB data analysts 
persistently claim that the inflation theory is in spectacular agreement with observations and the 
CMB data are perfectly well consistent with r = (the published confldence level contours always 
include r = and are typically centered at that point). 

Our analysis in this paper, based on general relativity and quantum mechanics, is aimed at 
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showing that there is evidence of signatures of rehc gravitational waves in the ah'eady available 
CMB data. We also make predictions for some future experiments and observations. 

The plan of the paper is as follows. In Sec. ^ we summarize the properties of a random 
background of relic gravitational waves that are needed for CMB calculations. The emphasis is on 
the gravitational wave (g.w.) mode functions, power spectra, and statistical relations. In Sec. II 111 
we discuss the general equations of radiative transfer and explain the existence of four invariants 
/, V, E, B that fully characterize the radiation field. We formulate the linearized equations in the 
presence of a single Fourier mode of gravitational waves. We prove that there exists a choice 
of variables that reduces the problem of temperature and polarization anisotropics to only two 
functions of time a(r/) and (3{ri). 

Section HVl is devoted to further analysis of the radiative transfer equations. The main result of 
this section is the reduction of coupled integro-differential equations to a single integral equation of 
Voltairre type. Essentially, the entire problem of the CMB polarization is reduced to a single inte- 
gral equation. This allows us to use simple analytical approximations and give transparent physical 
interpretation. In Sec. we generalize the analysis to a superposition of random Fourier modes 
with arbitrary wavevectors. We derive (and partially rederive the previously known) expressions for 
multipole coefficients af^ of the radiation field. We show that the statistical properties of the mul- 
tipole coefficients are fully determined by the statistical properties of the underlying gravitational 
perturbations. This section contains the expressions for general correlation and cross-correlation 
functions Cf-^' for invariants of the radiation field. 

We work out astrophysical applications in Sec. IVII and Sec. IVIII where we discuss the effects of 
recombination and reionization era, respectively. Although all our main conclusions are supported 
by exact numerical calculations, we show the origin of these conclusions and essentially derive 
them by developing and using semi-analytical approximations. The expected amplitudes, shapes, 
oscillatory features, etc. of all correlation functions as functions of ^ are under analytical control. 
The central point of this analysis is the TE correlation function. We show that, at lower multipoles 
^, the TE correlation function must be negative (anticorrelation), if it is induced by gravitational 
waves, and positive, if it is induced by density perturbations. We argue that this difference in sign 
of TE correlations can be a signature more valuable observationally than the presence or absence 
of the BB correlations. This is because the TE signal is about two orders of magnitude larger 
than the expected BB signal and is much easier to measure. We summarize the competing effects 
of density perturbations in Appendix IdI 

Theoretical findings are compared with observations in Sec. IVIIII In the context of relic gravi- 
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tational waves it is especially important that the WMAP team |9| stresses (even if for a different 
reason) the actual detection of the TE anticorrelation near £ ~ 30. We show that this is possible 
only in the presence of a significant amount of relic gravitational waves (within the framework of 
all other common assumptions). We analyze the CMB data and suggest models with significant 
amounts of gravitational waves that are consistent with the measured TT, TE and EE correlation 
functions. Our final conclusion is that there is evidence of the presence of relic gravitational waves 
in the already available CMB data, and further study of the TE correlation at lower ts has the 
potential of a firm positive answer. 

II. GENERAL PROPERTIES OF RELIC GRAVITATIONAL WAVES 

Here we summarize some properties of cosmological perturbations of quantum-mechanical ori- 
gin. We will need formulas from this summary for our further calculations. 

A. Basic definitions 

As usual (for more details, see, for example, Q|), we write the perturbed gravitational field of 
a flat FLRW universe in the form 

ds'^ = -c^dt^ + a^{t){6ij + hij)dx'dx^ = o?{'q) [-d-rf + {6ij + hij)dx'dx^] . (1) 

We denote the present moment of time hy rj = r]^ and define it by the observed quantities, for 
example, by today's value of the Hubble parameter Hq = H{r]ji). In addition, we take the present 
day value of the scale factor to be a(r?ij) = 21h, where Ih = c/Hq. 

The functions hij (7/,x) are expanded over spatial Fourier harmonics e^*"'^, where n is a dimen- 
sionless time independent wave-vector. The wavenumber n is n = {dijn'n^Y^'^ . The wavelength 
A, measured in units of laboratory standards, is related to n by A = 2-Ka/n. The waves whose 
wavelengths today are equal to today's Hubble radius carry the wavenumber tih = 47r. Shorter 
waves have larger n's and longer waves have smaller n's. 

The often used dimensional wavenumber defined by /c = 2ti /\[r]{i) in terms of today's wave- 
length X{r]ji), is related to n by a simple formula 

k = ^^n (1.66 X 10"'' h) Mpc"^ 

The expansion of the field hij {rj, x) over Fourier components n requires a specification of polar- 

s 

ization tensors Pij (n) (s = 1, 2). They have different forms depending on whether the functions 
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hij(rj,x) represent gravitational waves, rotational perturbations, or density perturbations. 

In the case of gravitational waves, two independent linear polarization states can be described 
by two real polarization tensors 

1 2 

Pij (n) = klj - rmmj, Pij (n) = kmj + rmlj, (2) 

where spatial vectors (1, m, n/n) are unit and mutually orthogonal vectors. The polarization tensors 
((^j) satisfy the conditions 

Pij5'^ = 0, kjn' = 0, iijp'^ = 25s's- (3) 

1 2 

The eigenvectors of Pij are li and mj, whereas the eigenvectors of Pij are li + nii and U — rrii. In 
both cases, the first eigenvector has the eigenvalue +1, whereas the second eigenvector has the 
eigenvalue —1. 

In terms of spherical coordinates 0, (/), we choose for (1, m, n/n) the right-handed triplet: 

1 = (cos cos cos6'sin</), — sin0), 

m=(— sin0, COS0, 0), (4) 

n/n = (sin cos 0, sin sin 0, cosO). 
The vector 1 points along a meridian in the direction of increasing 0, while the vector m points 
along a parallel in the direction of increasing (j). With this specification, polarization tensors @ 
will be called the '+' and 'x' polarisations. The eigenvectors of '+' polarization correspond to 
north-south and east- west directions, whereas the ' x ' polarization describes the directions rotated 
by 45°. 

With a fixed n, the choice of vectors 1, m given by Eq. (0]) is not unique. The vectors can be 
subject to continuous and discrete transformations. The continuous transformation is performed 
by a rotation of the pair 1, m in the plane orthogonal to n: 

1' = Icos^ -|- msin-;/', m' = — Isin^/; -|- mcosV', (5) 

where V' is an arbitrary angle. The discrete transformation is described by the flips of the 1, m 
vectors: 

1' = —1, m' = m or l' = 1, m' = — m. (6) 
When (jSJ is applied, polarization tensors transform as 



1 



1 2 



p',j (n) = I'il'j - m'im'j =P,j (n) cos 2^+ P,j (n) sin 2^^, 

2 12 ^ ' 

p'ij (n) = I'-m'j + m[l'- = - Pij (n) sin 2^--^ Pij (n) cos 2ip, 



and when (jEJ is applied they transform as 

1 1 2 2 

P'y (n) =Pi, (n), p'ij (n) = - Pij (n). (8) 

Later in this section and in Appendix ^ we are discussing the conditions under which the 
averaged observational properties of a random field hij{n, rf) are symmetric with respect to rotations 
around the axis n/n and with respect to mirror reflections of the axes. Formally, this is expressed as 
the requirement of symmetry of the g.w. field correlation functions with respect to transformations 
© and ©. 

In this paper we will also be dealing with density perturbations. In this case, the polarization 
tensors are 

These polarization tensors satisfy the last of the conditions ©■ 

In the rigorous quantum-mechanical version of the theory, the functions hij are quantum- 
mechanical operators. We write them in the form: 



+00 



kj (n) hn iv) e'"''' Cn +Pij (n)L (f?)e ^"''^ 



(10) 



where the annihilation and creation operators, Cn and c^, satisfy the relationships 

[c„, v]= W('^(n-n'), c„|0)=0. (11) 

The initial vacuum state |0) of perturbations is defined at some moment of time rjQ in the remote 
past, long before the onset of the process of superadiabatic amplification. This quantum state 
is maintained (in the Heisenberg picture) until now. For gravitational waves, the normalization 



constant is C = VlGvr/p;. 

The relationships 1)11(1 determine the expectation values and correlation functions of cosmolog- 
ical perturbations themselves, and also of the CMB's temperature and polarization anisotropies 
caused by these cosmological perturbations. In particular, the variance of metric perturbations is 
given by 



oo 

-'2 



f 2 r s dn 

(0| /..,(r/,x)/.*^(r/,x) |0) = ^ / | hn (r/)!^-. (12) 

i s=l,2 

The quantity 

h\n,r^) = f^n' ^ | L (^)P = ^ E I ^ (l^) 

^ s=l,2 s=l,2 



is called the metric power spectrum. Note that we have introduced 

h{n,ri) = -nhn{ri). (14) 
vr 

The quantity (|1H|) gives the mean-square value of the gravitational field perturbations in a logarith- 
mic interval of n. The spectrum of the root-mean-square (rms) amplitude /i(n, rj) is determined by 

___ s 

the square root of Eq. (|13|) . Having evolved the classical mode functions hn {'>]) up to some arbi- 
trary instant of time ij (for instance, today rj = rj^i) one can find the power spectrum h{n, rf) at that 
instant of time. For the today's spectrum in terms of frequency u measured in Hz, v = uHq/Att, 
we use the notation hymsi^)- 

In our further applications we will also need the power spectrum of the first time-derivative of 
metric perturbations: 



^ ' drj drj I ^ 2 y ^f^^ 



oo 

dn 







d h (n, T]) 



drj 



(15) 

n 



To simplify calculations, in the rest of the paper we will be using a 'classical' version of the 
theory, whereby the quantum-mechanical operators and are treated as classical random 
complex numbers Cn and c^. It is assumed that they satisfy the relationships analogous to Hll() : 

(Cn) = (c^,) = 0, (SnC„,) = (S„c„,) = dss'S^'^Hn - n'), {Lcn') = iKK') = 0, (16) 

where the averaging is performed over the ensemble of all possible realizations of the random field 

m- 

The relationships 1)16(1 are the only statistical assumptions that we make. They fully determine 
all the expectations values and correlation functions that we will calculate, both for cosmological 
perturbations and for the induced CMB fluctuations. For example, the metric power spectrum 
(|13|) follows now from the calculation: 

oo 

1 {h,,{v,x)h^^{v,x)) = I L m'^. (17) 



s=l,2 





s 

The quantities | hn {vn ^-re responsible for the magnitude of the mean-square fluctuations 
of the field in the corresponding polarization states s. In general, the assumption of statistical 
independence of two linear polarization components in one polarization basis is not equivalent to 
this assumption in another polarization basis. As we show in Appendix statistical properties 
are independent of the basis (i.e. independent of tp, Eq. ©), if the condition 

I hn (ry)P = I hn (r?)!' (18) 



is satisfied. As for the discrete transformations Q, they leave the g.w. field correlation functions 
unchanged. 

In our further discussion of the CMB polarization it will be convenient to use also the expansion 
of hij over circular, rather than linear, polarization states. In terms of definitions ((21), the left 
(s = L) and right (s = R) circular polarization states are described by the complex polarization 
tensors 

L * R R * L ^ ' 

Pij Piji Pij Piji 

satisfying the conditions (for s = L,R) 

kjS'^=0, Pijn' = 0, iijp'^* = 26s's. (20) 



A continuous transformation ((SJ brings the tensors H19|) to the form 



^ 4., e-^'^ P\,=i,e^'t (21) 



Functions transforming accor ding to the rule ((21]) are called the spin-weighted functions of spin +2 
and spin -2, respectively A discrete transformation ^ applied to (fT^ interchanges 

the left and right polarisation states: 

P'ij=Pij, P'ij=kj. (22) 

The assumption of statistical independence of two linear polarisation states is, in general, not 
equivalent to the assumption of statistical independence of two circular polarization states. More- 
over, it is shown in Appendix^that symmetry between left and right is violated, unless 

I hn (r?)!' = I hn (r?)|^ (23) 

However, if conditions (|18|) . (|23|) are satisfied, statistical properties of the random g.w. field remain 
unchanged under transformations from one basis to another, including the transitions between 
linear and circular polarizations. The summation over s in the power spectra such as ()17|) can be 

s 

replaced by the multiplicative factor 2. Moreover, the mode functions hn {f]) ^'^^ t'^o independent 
polarization states will be equal up to a constant complex factor e*°. This factor can be incor- 
porated in the redefinition of random coefficients without violating the statistical assumptions 
(|16j) . After this, the index s over the mode functions can be dropped: 

L (r?) = hniv)- (24) 



There is no special reason for the quantum-mechanical generating mechanism to prefer one 
polarization state over another. It is natural to assume that the conditions (|18|1 . (|23|) hold true for 
relic gravitational waves, but in general they could be violated. In calculations below we often use 
these equalities, but we do not enforce them without warning. 



B. Mode functions and power spectra 

s 

The perturbed Einstein equations give rise to the g.w. equation for the mode functions hn (v)- 
This equation can be transformed to the equation for a parametrically disturbed oscillator ^J: 



An + A„ 



a" 



0, (25) 



where An {f]) = ^(^) hn (??) and ' = d/d-q = {a/c)d/dt. (In this paper, we ignore anisotropic 
stresses. For the most recent account of this subject, which includes earlier references, see 
Clearly, the behaviour of the mode functions depend on the gravitational 'pumping' field a(ry), 
regardless of the physical nature of the matter sources driving the cosmological scale factor a(r/). 
Observational data about relic gravitational waves allow us to make direct inferences about H{rj) 
and a(?7) and it is only through extra assumptions that we can make inferences about such 
things as, say, the scalar field potential (if it is relevant at all). 



Previous analytical calculations (see, for example. 
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and references there) were based on 



models where a{'q) consists of pieces of power- law evolution 

«(r/)=^o|r/|i+^ (26) 

s f> 

where lo and j3 are constants. The functions a(ry), a'(r/), hn (??)i h'n (??) were continuously joined 
at the transition points between various power-law eras. 

It is often claimed in the literature that this method of joining the solutions is unreliable, unless 
the wavelength is "much longer than the time taken for the transition to take place" . Specifically, 
it is claimed that the joining procedure leads to huge errors in the g.w. power spectrum for short 
waves. It is important to show that these claims are incorrect. As an example, we will consider 
the transition between the radiation-dominated and matter-dominated eras. 

The exact scale factor, which accounts for the simultaneous presence of matter, oc a(r/)~^, 
and radiation, pj oc a(ry)~^, has the form 



ai,) = 2ln{\±^] V h + , (27) 



2 ~\~ Z(,q J \ 1 -|- Z(^q 
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where z^q is the redshift of the era of equahty of energy densities in matter and radiation Pm{zeq) 
For this model, the values of parameter tj at equality and today are given by the expressions 



^eq 



1 -\- Z^q 1 -j- Z^q 1 ~1~ Z^q 

The current observations favour the value Zeq^3x 10^ Q- 

The piece-wise approximation to the scale factor (|27|) has the form 

aiv) = r], r] < rjeq] a{r]) = 21h {v + Veqf , r/ > Veq, (28) 

Y -L + Zgq 



where rj^q and rm for this joined a(r/) are given by r^eg = l/2y/l + Z(,q and = 1 — 1/2-^/1 + Zeg. 
It can be seen from ()27() and H28() that the relative difference between the two scale factors is very 
small in the deep radiation-dominated era, r/ ^ l/^Ze<j) and in the deep matter-dominated era, 
r/ ^ 1/ ^Zeq- But the difference reaches about 25% at times near equality. 

The initial conditions for the g.w. equation (|25() are the same in the models ()27p . ()28() and 
they are determined by quantum-mechanical assumptions at the stage (which we call the i-stage) 
preceding the radiation epoch. The i-stage has finished, and the radiation-dominated stage has 
started, at some rji with a redshift Zi. The numerical value of Zi should be somewhere near 10^^ 
(see below). 

The initial conditions at the radiation-dominated stage can be specified at that early time r/j or, 
in practice, for numerical calculations, at much latter time, as long as the appropriate g.w. solution 
is taken as 

s s 

Pn {ff) = —2iB sinnr], p^ (?]) = —2inB cos nr], (29) 

where 



B = F{p) 



+ Zeq 
l + Zi 



and F[[3) is a slowly varying function of the constant parameter (3, |i<'(— 2)| = 2. Parameter (5 
describes the power-law evolution at the i-stage and determines the primordial spectral index n, 
n = 2/3 + 5. In particular, (3 = —2 corresponds to the flat (scale-invariant) primordial spectrum 
n = 1. 
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For numerical calculations we use the constant B in the form: 



B = 2 



47ry/l + z, 



1 + Zi 



■eq 



) 




(30) 



The wave-equation (|25() with the scale factor (|27() can not be solved in elementary functions. 
However, it can be solved numerically using the initial data (|29|1 . Concretely, we have imposed 
the initial data H29|) at r]r = 10^^, which corresponds to the redshift z ~ 3 x lO'^, and have chosen 
Zeq = 6 X 10^ for illustration. Numerical solutions for hn{ri) / hn{rjr) characterized by different 
wavenumbers n are shown by solid lines in Fig. ^ We should compare these solutions with the 
joined solutions found on the joined evolution ()28j) for the same wavenumbers. 

The piece-wise scale factor (|28() allows one to write down the piece-wise analytical solutions to 
the g.w. equation (^5]) : 



where 7/2 = nrjeq. In Fig. ^we show the joined solutions ()31() by dotted curves. 

It is clear from Fig. ^ that the g.w. solutions as functions of r] are pretty much similar in the 
two models. The solid and dotted curves are slightly different near equality (where the relative 
difference between the scale factors is noticeably large) and only for modes that entered the Hubble 
radius around or before equality. Moreover, the g.w. amplitudes of the modes that entered the 
Hubble radius before equality gradually equalize in course of later evolution. There is nothing like 
a huge underestimation or overestimation of the high-frequency g.w. power that was alleged to 
happen in the joined model. 

It is easy to understand these features. Let us start from waves that entered the Hubble radius 
well after equality, i.e. waves with wavenumbers n/riH ^ y^l + z^g. As long as these waves are 
outside the Hubble radius, their amplitudes remain constant and equal in the two models, despite 
the fact that the scale factors are somewhat different near equality. The waves start oscillating in 
the regime where the relative difference between (|27() and H28p is small, and therefore the mode 
functions, as functions of r], coincide. 




^Jn{r] + 'qeq) [A„ J3/2 ( n{ri + rjeq)) - iBnJ-3/2 ( ri{r] + 7?eg))] , V > Veq 



(31) 
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1.2 




time ri 

FIG. 1: The g.w. mode functions hn{vi) / hn{rir) in a matter-radiation universe. The solid curves are 
numerically calculated solutions on the scale factor (|27|l . while the dotted curves are analytical solutions on 
the scale factor (I28|l . 

The waves with wavenumbers njnu ^ ^\ + Zgg enter the Hubble radius well before the equal- 
ity. They oscillate in the WKB regime according to the law /i„(r/) oc e"^^"^ / a(r]) , having started 
with equal amplitudes in the two models. Near equality, the mode functions are different in the 
two models, as much as the scale factors are different. But the relative difference between ()27() and 
H28|) decreases with time, and therefore the mode functions in the two models gradually equalize. 
The amplitudes of these mode functions, as well as the scale factors H28() . are exactly equal 

today. The only difference between these high-frequency mode functions is in phase, that is, in 
different numbers of cycles that they experienced by today. This is because the moment of time 
defined as 'today' in the two models is given, in terms of the common parameter r], by slightly 
different values of r]ji. 

Finally, for intermediate wavenumbers n/nn ~ y^l + Zeq, the modes enter the Hubble radius 
when the scale factors differ the most. Therefore, they start oscillating with somewhat different 
amplitudes. The difference between these mode functions is noticeable by the redshift of decoupling 
Zdec (characterized by somewhat different values of rjdec), as shown in Fig.^ The difference survives 
until today, making the graph for the g.w. power spectrum (|1.S|) a little smoother (in comparison 
with that derived from the joined model) in the region of frequencies 10~^^ Hz that correspond to 
the era of equality. 
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frequency 



FIG. 2: The present-day spectra for hrms{v) and Vlgyj{v). The soHd hnes correspond to the primordial 
spectral index (3 = —1.9, i.e. n = 1.2, while the dashed lines are for j3 ~ —2, i.e. n = 1. 

Having demonstrated that the use of joined analytical solutions is well justified, and the previ- 
ously plotted graphs for hrmsi'^) ^-nd r2g^(i/) are essentially correct, we shall now exhibit the more 
accurate graphs based on numerical calculations with the initial conditions ()29() . H3U|) . We adopt 
Zeg = 3 X 10^ and Hq = 75^ /Mpc We adjust the remaining free constant Zi in such a way 
that the temperature correlation function £{i + l)Ci/2'iT at ^ = 2 is equal to 211 /zK^ This 
requires us to take Zi = 1.0 x 10^^ for f3 = —2 and Zi = 2.4 x 10^° for /3 = —1.9. 

The graphs for today's spectra, normalised as stated above, are shown in Fig. [2 In order to 
keep the figure uncluttered, the oscillations of hrms{t^) are shown only at low frequencies, while the 
oscillations of Qg^^i') are not shown at all. We have also indicated possible detection techniques 
in various frequency bands. 

The function Qgw{u) is the spectral value of the cosmological parameter [3]: 

\lgw{l'l,l'2) = — = — / Pgw[v) = / ilgw{v) , 

Pc Pc J Ul ^ Jui ^ 

that is, 

^gw{l^) = —Pgw{l^)- 
Pc 

Using the high-frequency definition of Pgw{^) (valid only for waves which are significantly shorter 
than Ih) we derive 



^gM = ^^hl^,[u)(^\\ (32) 



3 K'^hJ 
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It is this definition of 57^^(1/) that is used for drawing the curves in Fig. [21 

We have to warn the reader that a great deal of hterature on stochastic g.w. backgrounds uses 
the incorrect definition 



1 dpgu,{l>) 



Pc d Inu ' 

which suggests that the 17^^(1/) parameter is zero if the g.w. spectral energy density Pgw{^') is 
frequency-independent, regardless of the numerical value of Pgwi^^)- Then, from this incorrect 
definition, a formula similar to Eq. (|32() is often being derived by making further compensating 
errors. 

The higher- frequency part of hrms{v) is relevant for direct searches for relic gravitational waves, 
while the lower-frequency part is relevant to the CMB calculations that we turn to in the next 
section. The direct and indirect methods of detecting relic gravitational waves are considered in 



a large number of papers (see, for example, references 
though we disagree with some of them). 
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It is important to keep in mind that according to the infiationary theory the dashed lines for 
^gwi^^) and hrms{i^) in Fig. I^lshould be at a zero level, because they describe the g.w. background 
with a flat primordial spectrum /? = —2, n = 1. The 'consistency relation' of the inflationary 
theory demands r = 0, i.e. vanishingly small g.w. background, in the limit n = 1. 



III. EQUATIONS OF RADIATIVE TRANSFER 
A. Characterization of a radiation field 



The radiation field is usually characterized by the four Stokes parameters (/, Q, U, V), 23. l3C | . 
The parameter / is the total intensity of radiation, Q and U describe the magnitude and direction 
of linear polarization, and V is the circular polarization. The Stokes parameters can be viewed 
as functions of photons' coordinates and momenta Since photons propagate with the 

speed of light c, the momenta satisfy the condition PaP"' = 0. Thus, the Stokes parameters are 
functions of (t, x*, u, e*), where v is the photon's frequency, and e* is a unit vector in the direction 
of observation (opposite to the photon's propagation). 

In a given space-time point (t, x^), the Stokes parameters are functions of i', 0, (p, where 6, (p are 
coordinates on a unit sphere: 

da^ = gahdx^dJ' = dO'^ + sin^ Odc/)'^. (33) 
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The radial direction is the direction of observation. 

The Stokes parameters are components of the polarization tensor Pat jsol which can be written 

as 

1 / I + Q - (U-iV)sm9 \ 

\-{U + iV)sme {I-Q)sm^ej 

(We do not indicate the dependence of Stokes parameters on ly.) 

The symmetry of Pab with respect to rotations around the direction of observation requires the 
vanishing of linear polarization Q = 0, U = 0, but the circular polarization V can be present. The 
symmetry of Pab with respect to coordinate reflections in the observation plane requires = 0, 
but the linear polarization can be present. The first symmetry means that the readings of a linear 
polarimeter are the same when it is rotated in the observation plane. The second symmetry means 
that the readings of the left-handed and right-handed circular polarimeters are the same. (Compare 
with the discussion on gravitational waves in Sec. Ill Al ) 

Under arbitrary transformations of 9,(j), the components of Pab{0,(t>) transform as components 
of a tensor, but some quantities remain invariant. We want to build linear invariants from Pab 
and its derivatives, using the metric tensor Qabid, (f) and a completely antisymmetric pseudo-tensor 
e'^^(0,0). 



^ab 



- sin"^ 
sin-i e 



The first two invariants are easy to build: 

Iie,^)=g''\9,cP)Pabi9A), V{e,^) = ie''\e,^)Pabi9A)- (35) 



Then, it is convenient to single out the trace and antisymmetric parts of Pab, and introduce the 

ab 

1 



symmetric trace- free (STF) part Pa^^ ' 



Q -UsmO 



tdSTF 
^ab 



-Usm9 -Qsm'^9 



Clearly, the construction of other linear invariants requires the use of covariant derivatives of 
the tensor Paif^- There is no invariants that can be built from the first derivatives Pa^^ , so we 
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need to go to the second derivatives. One can check that there are only two linearly independent 
invariants that can be built from the second derivatives: 



:b:d 



(36) 



The quantities / and E are scalar s, while V and B are pseudoscalars. V and B change sign under 
flips of directions (coordinate transformations with negative determinants). This is also seen from 
the fact that their construction involves the pseudo-tensor e^fe. 

The invariant quantities (/, E, B), as functions oi 9,(j), can be expanded over ordinary spher- 
ical harmonics Yii„i{0 ^ (j)) ^ = 



1(0 A 
v{e,(i) 

E{9,4> 
B{0,(l) 



oo i 

f=0 m=-l 

oo I 



t=Q m=-e 
oo i 




(37a) 
(37b) 
(37c) 
(37d) 



The set of multipole coefficients {aj^, a^^, af^, af^) completely characterizes the radiation field. 
We will use these quantities in our further discussion. 

To make contact with previous work, we note that the multipole coefficients CLf^,af^ can also 
be expressed in terms of the tensor Pab itself, rather than its derivatives. This is possible because 
one can interchange the order of differentiation under the integrals that define a^, af^ in terms of 
the right hand side (r.h.s.) of Eq. This leads to the appearance of the spin-weighted spherical 



harmonics or tensor spherical harmonics 
written as 
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32l |. For example, the tensor Pab can be 



^ oo i 



le 



V 

ab'^im 



£=0 m 

oo 



+ 



i=2 m=-l 



where Y,9 



{£m)afe(^' and the "gradient" and "curl" tensor spherical harmonics 

forming a set of orthonormal functions for STF tensors The invariants E and B can also be 

n 

written in terms of the spin raising and lowering operators 9 and 3 |31l |: 



1 



E = -- [§2 (Q + iu) + g2 (Q _ iu)] , B = -[n^{Q + iu) - g2 (g - iu)] . 
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The quantities E and B are called the E (or "gradient") and B (or "curl") modes of polarization. 
The ^-dependent numerical coefficients in H37c|) and ()37d|) were introduced in order to make the 



definitions of this paper fully consistent with the previous literature 



31 
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B. Radiative transfer in a perturbed universe 



We need to work out the radiative transfer equation in a slightly perturbed universe described 
by Eq. The Thompson scattering of initially unpolarized light cannot gene rate circular polar- 



29 
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34J . we shall write the 



n 



1^ 

2/11/3 



I-Q 
-2U 



(38) 



ization, so we shall not consider the V Stokes parameter. Following 
radiative transfer equation in terms of a 3-component quantity (symbolic vector) n(x",p"). The 
components (ni,n2,ri3) are related to the Stokes parameters by 

(l + n\ 

1 C2 

where h is the Planck constant. The quantities fii,h2, (fii + n2 + n-3)/2 are the numbers of photons 
of frequency z/ coming from the direction z and passing through a slit oriented, respectively, in the 
directions x, y, and the bisecting direction between x and y. 

The equation of radiative transfer can be treated as a Boltzmann equation in a phase space. 
The general form of this equation is as follows 



35| 



as 



(39) 



where s is a parameter along the world-line of a photon, ^ is a total derivative along this world-line, 
and C is a collision term. We shall explain each term of this equation separately. 
The total derivative in Eq. (|39() reads: 



Dh 

ds 



dx" d dp"' d 
+ 



n. 



ds dx" ds dp" 

where dx" /ds and dp" /ds are determined by the light-like geodesic world-line. 



dx" 
ds 



P , 



dp" 
ds 



r" rJ^n'y n ort"n 
A P ) QafiP P 



0. 



(40) 



(41) 



Strictly speaking, the square bracket in Eq. (|4n|) should also include the additive matrix term R. 
This term is responsible for the rotation of polarization axes that may take place in course of 



parallel transport along the photon's geodesic line 



3q |. In the perturbation theory that we are 
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working with, this matrix does not enter the equations in the zeroth and first order approximations 
[3^, and therefore we neglect R. 

In our problem, the collision term C describes the Thompson scattering of light on free (not 
combined in atoms) electrons 29(. We assume that the electrons are at rest with respect to one 
of synchronous coordinate systems (^. We work with this coordinate system, so that it is not 
only synchronous but also 'comoving' with the electrons. (This choice is always possible when 
the functions hij in ^ are gravitational waves. Certain complications in the case of density 
perturbations will be considered later, Appendix IdI) Thus, the collision term C is given by the 
expression 



C [n] = -arNeix'^) (^^^ h{t, x\ v,e,<P)-^ j dn'P{e, 0; 9', (l)')h{t, x\ v, e\ 0') 



(42) 



where ax = 6.65 • 10^^^ cm^ is the Thompson cross section, N^, is the density of free electrons, and 
F{6,(j);6' is the Chandrasekhar scattering matrix. (The explicit form of the scattering matrix 
is discussed in Appendix lO) The factor cdt/ds arises because of our use of the element ds, instead 
of cdt, in the left hand side (l.h.s.) of Eq. H39() . In accord with the meaning of the scattering term, 
the quantity aTNe{x°'){cdt / ds) is the averaged number of electrons that could participate in the 
scattering process when the photon traversed the element ds along its path. 

Let us now write down the equations of radiative transfer in the presence of the gravitational 
field (^. First, we write down the equations for the light-like geodesic line a;"(s) = (77(5), a;*(s)): 



r, dn u dx* v ,• du 

p^ = -r = —, P =-r = — e ' -r 

ds ca ds ca ds 



{5ij + hij) e'e^ = 1. 



1 da 1 A dhij 

1 — e e — 

adrj 2 di] 



We do not need the expression for de^ /ds, because it is a first-order (in terms of metric perturba- 
tions) quantity, and this quantity enters the equations of radiative transfer only in products with 
other first-order terms. We neglect such second-order corrections. 
Second, we write for the 'vector' n: 

n = n(o)+n«, (44) 

where 

n(0) is 

the zeroth order solution, and fi^^) is the first order correction arising because of 
the presence of metric perturbations. We shall now formulate the equation for n^^^, taking into 
account the zero-order solution to Eq. (jSH). 
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In the zero-order approximation we assume that hij = and that the radiation field is fully 
homogeneous, isotropic, and unpolarized. Therefore, 



n 



W =nW(r?,i/)u, 



(45) 



where 

u = 1 

vv 

Since the scattering matrix P does not couple to the radiation field if it has no quadrupole 
anisotropy, the collision term H42|) vanishes, C [n'-*'-*] = 0, and the equation for n^^\rj,i') reads 

9n(°) u da an(0) _ ^ 
dr] a dr] dv 

The general solution to this equation is n^'^'> = uq {1/0(7])), which makes it convenient to use a new 
variable 

D = va{rj). 

In the zero-order approximation, v = const /a{rj) and therefore i> is a constant along the light ray. 

We are now in a position to write down the first-order approximation to the Boltzmann equation 
((Snj- We take {r],x\i>,e^) as independent variables, i.e. n^^^ = no{h'), n.^^^ = n^-^^ (r/, x*, i>, e*), and 
use the identity 

dp" d di> d ^de^ d 
ds dp" ds du ds de^ 

in the first-order approximation to ()4U() . Taking also into account the geodesic equation H43|l we 

arrive at the equation 



h e : ve. e' 

dr] dx^ 2 df] dv 



^ = C 

ds 



fid) 



Introducing new notations g(ry) = aTOi{r])N e,(r]) and f(y) = (?ln no/(? In t' (the astrophysical 
meaning and numerical values of the functions qij]) and f(y) are discussed in Appendix El) we 
write down the final form of the transfer equation: 



d ^ ^^^^ _^ g« ^ 
dr] dx^ 
f{i')no{i>) j ,dhij , 

2 d^ 



u + g(r?)^ y P(e^e'^)n(i)(r?,x\i>,e'-'). 



(46) 
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It is seen from Eq. (|46() that the 'source' for the generation of n*-^-* consists of two terms on 
the r.h.s. of this equation. First, it is the gravitational field perturbation hij, participating in the 
combination e^e^dhij/drj. It directly generates a structure proportional to u, i.e. a variation in 
the / Stokes parameter and a temperature anisotropy. In this process, a quadrupole component of 
the temperature anisotropy necessarily arises due to the presence of the term e^d/dx^, even if the 
above-mentioned combination itself does not have angular dependence. The second term on the 
r.h.s. of Eq. H46|) generates polarization, i.e. a structure different from u. This happens because 
of the mixing of different components of h^^\ including those proportional to u, in the product 
term Pn^^^. In other words, polarization is generated by the scattering of anisotropic radiation 
field. Clearly, polarization is generated only in the intervals of time when q{r]) ^ 0, i.e. when free 
electrons are available for the Thompson scattering (see, for example, js^). 



C. The radiative transfer equations for a single gravitational wave 

We work with a random gravitational field hij expanded over spatial Fourier components (|1U() . 
It is convenient to make similar expansion for the quantities n^^\r],x^,i',e'^). Since Eq. (|46() is 
linear, the Fourier components of n^^^ inherit the same random coefficients Cn that enter Eq. ()1U() : 

+ CX) 



(27r)3/2 



/2n 



s=1.2 



(47) 



Equation (|46jl for a particular Fourier component takes the form: 



dr] 



+ Qiv) + ^eVj 



= l^^e^e^ P, („)%Mu + Mfdn' P(e^ e'^)n« (., ^, e'^- (48) 
2 dr] 47r J 

To simplify technical details, we start with a single gravitational wave propagating exactly in 
the direction of z in terms of definitions 0, @, (|19jl . The coordinate system, associated with 
the wave, is specified by = 0, = 0. This simplifies the polarization tensors H19|) and makes 
them constant matrices. At the same time, the observation direction is arbitrary and is defined 
by e* = (sin0cos(/), sin sine/;, cos 9). We consider circularly polarized states with s = 1 = L, 
s = 2 = R. Then, we find 



e'ei Pij (n) = (1 - ]i^)e 



(49) 
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where ^ = cos 9 and the it signs correspond to s = L,R, respectively. This simpUfication of the 
angular dependence is possible only for one Fourier component, but not for all of them together. 
We shall still need the results for a wave propagating in an arbitrary direction. The necessary 
generalization will be done in Sec. IV Bl 

The ziz2(j) angular dependence of the source term in Eq. (|48|) generates the ±2(1) angular depen- 
dence in the solution 23, js^. We show in Appendix [0 that the terms in h^]i{r],i>, with any 
other ^-dependence satisfy homogeneous differential equations, and therefore they vanish if they 
were not present initially (which we always assume). Similarly, the u dependence of the solution 
can be factored out. Finally, we show in Appendix O that one linear combination of the three 
components of fin ^ always satisfies a homogeneous equation and therefore vanishes at zero initial 
data. Thus, the problem of solving Eq. (|15|l reduces to the problem of finding two functions of the 
arguments r], fj,. 

Explicitly, we can now write 

(l\ ( (1 + /^') ^ 



1 



(50) 



Clearly, function a is responsible for temperature anisotropy (I Stokes parameter), while function 
f3 is responsible for polarization {Q and U Stokes parameters). 



Temporarily dropping out the labels n, s and introducing the auxi^ 



a{r],fi) + P{ri,ij), we get from Eqs. (|i8|) . (|50|) a pair of coupled equations 34l |: 



iary function £,{rj,fi) 



dr] 



where 



1 



dr] 



(l + ^'Y/3(7?,^')-2(l-^")'^(^'/^') 



(51) 
(52) 

(53) 



IV. RADIATIVE TRANSFER EQUATIONS AS A SINGLE INTEGRAL EQUATION 

In some previous studies [^, 0, 0|, equations (|^ . are being solved by first expanding 
the /^-dependence of /3 and ^ in terms of Legendre polynomials. This generates an infinite series of 
coupled ordinary differential equations. Then, the series is being truncated at some order. 
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We go by a different road. We demonstrate that the problem can be reduced to a single 
mathematically consistent integral equation. There are technical and interpretational advantages 
in this approach. The integral equation enables us to derive physically transparent analytical 
solutions and make reliable estimates of the generated polarization. The numerical implementation 
of the integral equation considerably saves time and allows simple control of accuracy. 



A. Derivation of the integral equation 



In order to show that the solutions of Eqs. (|51jl . (|52|) for a{r],fj.) and (3(ri, fi) are completely 



determined by a single integral equation, we first introduce new quantities 



4C 



lb 



(54) 



H{r]) = e 



Solutions to Eqs. (j^T|) . can be written as 



-T(n) dhji]) 
dt] 



(55) 



(56) 



(57) 



Expression H56|) is a formal solution to Eq. H51|) in the sense that P{rj,fj,) is expressed in terms of 
$(ry) which itself depends on (3{r],fi) (see (|5:-{|) and (jSU)- 

We now put (|56|) and (|57j) into Eq. (|53|) to get a new formulation for X(r/) : 



1 r, 

Z{rj) = e^^'') j j dfidr]' 

-1 

Using the kernels K±{rj — rj'), 



(58) 



K±irj-V') = I rfMli^') 



(59) 



Eq. (|5S|) can be rewritten as 



X(r]) = e^('') j dri 




K^(j] - ) - -K^{r] - r]')H{r]') 



(60) 
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Multiplying both sides of this equality by (3/16)g(ry)e '^^'^^ and recalling the definition ((511) we 
arrive at a closed form equation for ^{r]): 

■n 

HV) = ^q{ri) I drj'^rj')K+{r, - r,') + F{rj), (61) 


where F(r]) is the known gravitational-field term given by the metric perturbations, 

F{ri) = -^Qiv) I dri'H{rj')K^{r] - r?'), (62) 



The derived equation (|61() for ^{rj) is the integral equation of Voltairre type. As soon as <&(??) is 
found from this equation, we can find /3(ry, /i) from Eq. (|56|) . Then, Eqs. (|57|) and (|56|) completely 
determine all the components of n^^^ according to Eq. (j50() . 

Clearly, we are mainly interested in temperature and polarization anisotropics seen at the 
present time rj = rjj^. Introducing ( = n{r]R — rj) and restoring the indices n and s, we obtain the 
present-day values of a and j3: 

VR. 

an,s{lj) = an,s{VR, fJ-) = y [Hn,s{v) - ^n,si'n)^ c"*'"^, (63a) 



Vr 

PnAf^) = PnAm,f^) = I dr, $„,,(r/)e-^'^^. (63b) 



The integrals can safely be taken from 77 = as the optical depth r quickly becomes very large in 
the early Universe, and the 'source' functions Hn^siv) ^^^^1 <I>„^s(r/) quickly vanish there. We will 
work with expressions (|63|) in our further calculations. 



B. Analytical solution to the integral equation 

The integral equation (j6H) can be solved analytically in the form of a series expansion. Although 
our graphs and physical conclusions in this paper are based on the exact numerical solution to Eq. 
(|61|) . it is important to have a simple analytical approximation to the exact numerical solution. 
We will show below why the infinite series can be accurately approximated by its first term and 
how this simplification helps in physical understanding of the derived numerical results. 

We start with the transformation of kernels ()59() of the integral equation ()61|) . Using the identity 
Ij^k^ix/i _ (^(}_li(lx)^e^^^ , the kernels can be written as 

+1 2 
KAV-V) = jd,{l±,'fe-^^^-^') =2(lT^) 
-1 
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where x = n{r] — rj'). Now, taking into account the expansion 



smx 



X 



E 

m=0 



7X , 



(2m + 1)!' 



the r.h.s. of Eq. (jblj) can be presented in the form of a series: 



^r^) = -q{r]) n^" / dr]' {rj - r^'f'' A+(m)$(r/') - \.{m)H{r]' 



where 



A+(m) 



(-ir 



(2m + 1)! 



m=0 



1-4 



(64) 



(m + 2) 



(2m + 3)(2m + 5) 



A-_(m) 



-1) 



(2m + l)!(2m + 3)(2m + 5) ' 



Since the r.h.s. of Eq. (|64|) is a series in even powers of the wavenumber n, the l.h.s. of the 
same equation can also be expanded in powers of v?'^: 



,2m 



(65) 



m=0 



Using expansion H65|) in both sides of Eq. H64|) we transform this equation to 
J]$(-)(ry)n2™= -g(r/) 



m=0 



^ ^ A+(m)n 

m=0 j=0 



2{m+j) 



j dri'^^^\ri'){'q-ri'] 


^ A_ (m)n2'" y dr]'H{ri'){r] - rjf' 



2m 



(66) 



The left side and the second term in the right side of Eq. (|66() are series in n?^, but the first sum 
on the r.h.s. of Eq. (|66|) is still a mixture of different powers. This sum can be rearranged to be 
manifestly a series in n^"*: 

V 



^ ^ A+(m)n2('^+j) / d7?'$(j)(r?')(r/ - r/'; 



2m 



m=0 j=0 



m=0 



2m 



f;A+(A:) '=)(r?')(7?-V; 
fc=o ;( 



2A: 



(67) 



According to Eq. ()66() we have to make equal the coefficients of terms with the same power n^™' 
in both sides of the equation. This produces a set of integral equations 



(68) 
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where functions S^'^\rj) depend only on the known function H{r]) and functions ^^"^ ''^v) Pi'^- 
sumed to be found from equations of previous orders: 

S(-"-^)(r,) = -^A_(m) j dT]' H{r]'){r]-i]'f"' 


V 

+^f;^+(^) / ^V' ^^"'-'Hv'Xv-vT- (69) 
k=i i 

The important advantage of the performed expansion in powers of n^™ is that Eq. (|68|) of any 
order m is now a self-contained analytically solvable integral equation. 

Exact solution to the integral equation H68|) is given by the formula 

v 

J df]' 


We can further simplify this formula. Taking an //-derivative of expression (|69|1 we find 



J\2m-1 



^ = -Y^H{rj)dom - uj-{m) J drj H{t] ){v-V J 





k=l 

where uj±{m) = 3mX±{k). Substituting this expression into Eq. H70|). we arrive at the final result 

V 







(71) 



+ 5]u;+(A;)cI>(— ^■)(^')^(fe)(^,r/' 

k=l 

where 

*W(^,^0 = / dv" e-i^-^^"'^'Hv" - VT'-'- (72) 
n' 

Functions ^ (^j^^[r],'q') depend only on the ionization history of the background cosmological 
model described by q{ri). These functions can be computed in advance. Complete determination 
of the function <I>(r?), Eq. requires only one integration by r/ at each level m in Eq. (|71|) . 

starting from m = 0. The zero-order term (^^^\'q) does not depend on functions '^(k) and is 
determined exclusively by H{i]), Eq. ()55p. The zero-order term can be presented as 

V 

(r?) = -^^9{V) I drj'^!^e-^o^(^'^'\ (73) 
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It is crucial to remember that the function $(ry), Eq. (|54() . always contains the narrow visibility 
function g(r]) (see Appendix IH)) . In particular, function ^^^\ri) is nonzero only for rj within 
the width of g{r]), and is proportional to this width. In the era of decoupling, we denote the 
characteristic width of g{r]) by Ary^ec- With Aijdec we associate the characteristic wavenumber n*: 

2tt 
= -i: • 

Numerically, Aj^^gc ~ 3 x 10~^ and n,,, ~ 2 x 10^. In what follows, we will be interested in CMB 
multipoles i < 10'^. They are mostly generated by perturbations with wavenumbers n < 10^. 
Therefore, for wavenumbers of interest, we regard n/n^: as a small parameter. 

We shall now show that <I>(*^)(r/) is the dominant term of the series ()65() . The next term, 
<I>(^)(?7)n^, is at least a factor (n/n^,)^ smaller than <I>(*^)(r/), and so on. The explicit expression for 
^^^^ is as follows 



Effectively, function (^^^\ri), in comparison with ^>W(7?), contains an extra factor ^'(^^^(r], r]'). Tak- 
ing into account the fact that the functions g{ri) and e localized in the interval of 
arguments not larger than Arjdec, this factor evaluates to a number not larger than (Ar/f^ec)^- 
Therefore, the m = 1 term in Eq. H65|) is at least a factor (n/n^)'^ smaller than the m = term. 

These analytical evaluations are confirmed by numerical analysis as shown in Fig. |5J The solid 
line shows the exact numerical solution found from Eq. ()61() for hn{ri) and q{rj) described in Sec. Ill Bl 
and Appendix^ respectively. The dashed line is plotted according to formula ((75|) . with the same 
/i„(r/) and q{r]). It is seen from Fig. EJthat the zero-order term <I>('^)(r/) is a good approximation. 
The deviations are significant, they reach (20-25)%, only for the largest wavenumbers n in the 
domain of our interest. 



■(1 



dr]' 



(l)e^(V)$(0)(^/) 



■^-('''^')*(i)(r?,r?' 



V. MULTIPOLE EXPANSION AND POWER SPECTRA OF THE RADIATION FIELD 
A. Multipole coefficients 

Having found ^n,si'n) for a single gravitational wave specified by Eq. H49|l one can find a and 
P functions according to Eqs. (jlvEj) . (|63b|) . Then, using Eqs. (jMl), (El), ^ and (jHEI), one 

can find the multipole coefficients {X = I, E, B) participating in the decompositions pT]) . 
Although this route has been traversed before |3l|, we have made independent calculations in a 
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FIG. 3: Function ^n{v) for different values oin. The solid line is exact numerical solution to Eq. The 
dashed line is the zero order approximation (|73|l . The dotted line shows the approximation (|91(l (see below). 
The g.w. mode functions are normalized such that hn{rir) = 1. 

more general arrangement. Formulas derived in this subsection are effectively a confirmation of 
the correctness of calculations in Ref. ,31] • 

First, we integrate over photon frequencies (fot the definition of 7 see Eq. dHll)) and arrive at 
the following expressions 



En,s (/", </>) 
Bn,s (/", <t>) 



-7 



(1 - f^') ((1 + /^') 1^ + 8^1^ + 12) f3n,s e 



(74a) 
(74b) 

(74c) 



where, as before, the upper and lower signs correspond to s = 1 = L and s = 2 = R, respectively. 
The ziz2(j) dependence in Eqs. (|74j) implies that only the m = ±2 multipoles are nonzero. 

Then, we integrate Eqs. (fTi]) over angular variables in order to find according to Eq. pTj) . 
Using the notations T£{(), Ei((), Bi{() for the functions arising in course of calculations (and 
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called multipole projection functions) 

m) = ^ 

EdO = 
BeiC) = 



' + 2)! MC) 
'-2)\ C ' 
l{l - 1 



c 



MO+je-iiO 



and replacing a(/x) and by their expressions (|63|) . we finally arrive at 

aj^{n,s) = i-iY~'^ (^52,mSi,s + S-.2,mS2,s^ af {n, s), 

af^{n,s) = {-iY~^ (^S2,mSi,s - 6-2,mh,s^ af {n, s), 



where 



Ti 

at 



{n,s) = 7 V47r(2£ + 1) / dr, (Hr^ - ^-n.s (??)) T>(C), 





VR 



\n,s) = ^^^T^(2i + 1) / $„,,(r/)^,(C), 





Vr 



\n,s) = 7 V4vr(2£ + 1) ^ dr] cI>„,,(r/)i?^(C). 



(75a) 
(75b) 
(75c) 

(76a) 
(76b) 
(76c) 

(77a) 
(77b) 
(77c) 



B. Superposition of gravitational waves with arbitrary wavevectors 

It is important to remember that the result (|76|) is valid only for a special wave, with the 
wavevector n oriented exactly along the coordinate axis z. Since the perturbed gravitational 
field is a random collection of waves with all possible wavevectors n, and we are interested in 
their summarized effect as seen in some fixed observational direction ^, </), we have to find the 
generalization of Eq. H76|) to an arbitrary wave, and then to sum them up. 

To find the effect of an arbitrary wave, there is no need to do new calculations. It is convenient 
to treat calculations in Sec. IV Al as done in a (primed) coordinate system specially adjusted to a 
given wave in such a manner that the wave propagates along n' = (0,0, n). The observational 
direction e* is characterized by 0',(f)'. The quantities X = I,E,B calculated in Sec. IV Al are 
functions of 9',(l)' expanded over Yim{d' , cp'), 



(78) 



=0 m= 
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where af^ is a set of coefficients (compare with Eq. p7])l: a^, [(£ + 2)!/(£ — 2)!]^''^ a^, 
[(£ + 2)!/(^-2)!]i/2af„. 

Now, imagine that this special (primed) coordinate system is rotated with respect to the ob- 
server's (unprimed) coordinate system by some Euler angles 



a = (t>n, P = 0n, 7 = 



(see, for example 



4l|). The same observational direction e* is now characterized by 9,(1), and the 



same wavevector n' is now characterized by the unit vector 

n = n/n = (sin 9^ cos sin^n sin (p^, cos ^n)- 

Obviously, the already calculated numerical values of the invariants X{9',(f)') do not depend on 
the rotation of the coordinate system. Being expressed in terms oi 9,(j), the invariants describe the 
effect produced by a wave with a given (arbitrary) unit wavevector n, as seen in the direction 9, cj). 

The transformation between coordinate systems {9' , cf)') and {9, (p) is accompanied by the trans- 
formation of spherical harmonics, 

e 

Yt^(9',4l)= Di,^^{h)Y,rn'{9,<t>), 

m'=-l 

where 

^m',m (n) = r)ln<,m {<Pn,9n,0) 

are the Wigner symbols [4^. Later, we will need their orthogonality relationship 

J dnDl^pih)Dil^in) = J^5u'Smm' , (79) 

where = sin 9nd9nd<pn- 

We can now rewrite Eq. ((TSj) in terms of 9, (p, and thus find the contribution of a single arbitrary 
Fourier component. 



00 e / e \ 
e=o m=-e \m'=-e J 



(80) 



The superposition of all Fourier components of the perturbed gravitational field gives, at the 
observer's position x = and at t] = rjji (see (|47p ) . the final result: 



,~oo 

X{9,<P) = / ^ E [XnAO,cP) +X*^A9,<P) t] . (81) 



-00 



s=l,2 
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From this expression, combined with Eq. (|8U() . one can read off the random multipole coefficients 
afj^ that participate in the expansions (|37|): 




(82) 

Since X is a real field, the multipole coefficients obey the reality conditions 

as is seen directly from H82() . The gravitational- wave nature of metric perturbations is encoded in 
concrete values of the coefficients (f77)) . (|^ . But in all other aspects the argumentation presented 
here is general. 



C. Angular power spectra for temperature and polarization anisotropies 

It follows from Eq. (|82|) that the statistical properties of the multipole coefficients are fully 
determined by the statistical properties of the gravitational field perturbations represented by the 
random coefficients c^. A particular realization of is responsible for the particular realization 
of af^ actually observed in the sky. Having derived the distribution function for from some 
fundamental considerations (for example, from the assumption of the initial quantum-mechanical 
vacuum state of perturbations) we could estimate the probability of the observed set af^ within the 
ensemble of all possible sets. We could also evaluate the inevitable uncertainty in the observational 
determination of the parameters of the underlying random process. This uncertainty is associated 
with the inherent absence of ergodicity of any random process on a 2-sphere (i.e. sky) I4J]. In this 
paper, however, we adopt a minimalistic approach; we postulate only the relationships H16|) and 
calculate only the quadratic correlation functions for a^. 

Clearly, the mean values of the multipole coefficients are zeros, 

To calculate the variances and cross-correlation functions, we have to form the products af^a^^, 
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and then take their statistical averages. First, we find 

f n^dndQ 



/ X* X' \ ^ f n dndVt y—^ ^—^ ^—^ 

(a^™«^'m') = 7^ y 2n 2^ 2^ 2^ 



s=\,2mY=-lra\=-l 

'4m^ {n, s)a^^, (n, s)D^^^^ (n)I)^/„/ (n) 

{ri, s)af^,^ (n, s)D'_^^„,^ {^)D% ^„,.^ (fi)] . (84) 

We now take into account the fact (compare with Eq. ()76p) that 

«ftni iji, s)a^^, (n, s) oc 5^^^/ . 

This property allows us to get rid of summation over m'l in Eq. H84j) . Then, we perform integration 
over dJl and use the orthogonality relationships (|79() . We finally arrive at 

{af^a^^,) = See' 5mm', (85) 

where 



Cf^' = 4^2(2^ + 1) J '^dnY^ [afm{n,s)afm{n,s) + afm{n,s)afm{n,s)j . (86) 



Other quadratic averages, such as {O'fm^^m')^ i^fm'^em')^ follow from and the reality condition 
(jHl. 

The angular correlation and cross-correlation functions of the fields /, E, B are directly express- 
ible in terms of Eq. (|8H|) . For example, 

(/(^i,<Ai)/(^2,<A2)) = r(<5) = ^_^cr^'£(cos5), 

£=0 

where 5 is the angular separation between the directions (0i,(/)i) and {02,4>2) on the sky. If the 
actually measured values of af^ represent a particular realization of a gaussian random process, 
quantities Cf^', constructed from the measured af^ according to the r.h.s. of Eq. (|HH|l. are the 
best unbiased estimates 

One can note that the final result (|85() , ()86() is the integral of individual contributions (|76() from 
single gravitational waves given in a special frame discussed in Sec. IV A[ However, one cannot 
jump directly from (|76|1 to (|86|) (which is often done in the literature). In general, Eq. (|86|) does 
not follow from Eq. ()76|). By calculations given in this subsection we have rigorously shown that 
Eq. (|H?t|) is justified only if special statistical assumptions (|TH|) are adopted. 

One can also note that the correlation functions containing the label B once, i.e. CP and Cf ^, 
vanish if the extra assumptions ((TH|l. are made. Indeed, under these assumptions one can use 
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FIG. 4: The 'source' functions Hn.s{f]) ^^^d ^n.s{v) — 100) of temperature and polarization anisotropics 
(the normalization is chosen such that hn{rjr) = 1). 

one and the same mode function for both polarization states s, Eq. ()24() . Then, the coefficients 
a^(n, s), Eq. H76c|) . differ essentially only in sign for two different s, i.e. o£(n,-L) = —af_2{n,R). 
Therefore, their contributions will cancel out in expressions (jHEI) for Cf^ and Cf^. This statement 
is in agreement with Ref. \^\ . 

Without having access to a 'theory of everything' which could predict one unique distribution 
of the CMB radiation field over the sky, we have to rely on the calculated statistical averages (j86|l . 
We can also hope that our universe is a 'typical' one, so that the observed values of the correlation 
functions should not deviate too much from the statistical mean values. 

VI. EFFECTS OF RECOMBINATION ERA 

All our final graphs and physical conclusions in this paper are based on exact formulas and 
numerical calculations, starting from numerical representation of the key functions H(r]) and 
Eqs. (|55|) . (|6H) . However, we derive and explain all our results by developing manageable and 
accurate analytical approximations. At every level of calculations we compare exact numerical 
results with analytical ones. 

A. Temperature anisotropy angular power spectrum 

The temperature anisotropy power spectrum Cj^ is determined by the multipole coefficients 
aj{n,s), Eqs. (|86|) . (|76aj) . (|77aj) . The typical graphs for the functions Hn^siv) ^^'^ ^n,s{'n) 
shown in Fig. [l] 

Since the visibility function g{r]) is a narrow function, a convenient analytical approximation 
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is the hmit of an instanteneous recombination. The function e~'^ is replaced by a step function 
changing from to 1 at = Tj^ec, e""^ = h{r] — Tj^ec), and the function g{7]) is replaced by a delta- 
function, g{r]) = 6{r] — "qdec)- In this limit, the contribution to aj{n,s) from the scattering term 
^n^siv) is proportional to Arjdec- It can be neglected in comparison with the contribution from the 
gravitational term Hn^siv)- The ratio of these contributions is of the order of nAr]rec, and it tends 
to zero in the limit of instanteneous recombination, Aij^ec 0. 
Neglecting the scattering term, we write 

Vr s 

aj{n,s) = 7V4vr(2^ + l) / dv ^m). (87) 

Vdec 

This integral can be taken by parts, 

VR ^ 
hn (.'ndec)Te{Cdec) + J dj] n hn (^)~^^ 

Vdec 

The remaining integral contains oscillating functions and its value is smaller, for sufficiently large 
n's, than the value of the integrated term. This is illustrated in Fig. O Therefore, we have 



aj(n,s) =7V47r(2£ + l) 



aj{n, S) = -7^471(2^ + 1) hn (%ec)7>(Cdec). (88) 

Finally, we put Eq. (|S5|) into Eq. and take into account the definition of the metric power 
spectrum (|13p. Then, we get 

Cr = W I ^h\n, Vdec)Ti{Cdec). (89) 

The projection factor Tf{Cdec) is given by Eq. ()75a() . Since the spherical Bessel functions 
reach maximum when the argument and the index are approximately equal, Cdec ~ ^, a partic- 
ular wavenumber n is predominantly projected onto the multipole d k. n: 

^ ~ Cdec = ^{'qR - Vdec) ~ n. (90) 

This can also be seen in Fig. [51 Thus, the oscillatory features of the metric power spectrum 
h'^{n,r]dec) in the n-space are fully responsible for the oscillatory features of the angular power 
spectrum £{£ + \)Cj'^ in the ^-space (We use this opportunity to correct a misprint in Fig. 2 
of Ref. 1]| : the plotted lines are functions C;, not l{l + l)Ci. For the early graphs of Ci see Ref. 

The g.w. metric power spectrum h'^{n,rjdec) for the case /3 = —2, and the function £{£ + 1)C'['^ 
caused by this spectrum, are shown in Fig. [Tfa) and Fig. [7|^b). The normalization of the metric 
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FIG. 5: The contributions to the power spectrum t{l + l)Cj"^ from an individual mode n. The solid line 
shows the exact result calculated according to (I87II . while the dashed line shows the approximation H88|l . 
The normalization has been chosen such that hn(r\r) = 1. 

power spectrum is such that the function 1(1 + l)Cj'^ at ^ = 2 is equal to 1326 /iK^ The 
interval ^ < 90 is generated by waves with n < 90. These waves did not enter the Hubble radius by 
the time r/^ec- Their amplitudes are approximately equal for all n's in this interval (compare with 
Fig. The gradual decrease of the angular power spectrum at larger ^'s is the reflection of the 
gradual decrease of power in shorter gravitational waves whose amplitudes have been adiabatically 
decreasing since the earlier times when the waves entered the Hubble radius. 

B. Polarization anisotropy angular power spectrum 

The decisive function for polarization calculations is ^n,s(j])- We have approximated this func- 
tion by ^n'ls, Eq. (|73() . and compared it with exact result in Fig. |5J For qualitative derivations it 
is useful to make further simplifications. 

Since g(r\) is a narrow function, the integral in Eq. (|73|) is effective only within a narrow interval 

s 

Ar/rfec- Assuming that the function d hn {v)/dri does not vary significantly within this interval, we 
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can take this function from under the integral 



(91) 



Clearly, the assumption that the function d hn (^)/'^?? is almost constant within the window 
^Vdec gets violated for sufficiently short waves. With 1)91^ . we expect degradation of accuracy for 
wavenumbers n approaching n*. This is illustrated by a dotted line in Fig. |5J Nevertheless, the 
approximation (|91() is robust for n < n*, and it reveals the importance of first derivatives of metric 
perturbations for evaluation of the CMB polarization. 

We now introduce the symbol P to denote either E ox B components of polarization. In terms 
of multipoles (n, s), the angular power spectrum is given by Eq. l|86|) . Putting (|9T|l into Eqs. 
(|77b|) . H77cl) and denoting by ^^(C) the respective projection functions, we get 



af in, s) 



7 V47r(2£ + 1) j c 



dr] 



I dhn jr]) 
10 dr] 



PiiO 







(92) 



Again referring to the peaked character of g{rf) and assuming that the combination {d hn 
{f]) / dr})Pi[C,) does not change significantly within the window A-qdea we take this combination 
from under the integral, 

/ s 

I dhn (ll) 



af(n,s) « 7V47r(2£ + l) D{ 



n] 



10 dr] 



-m) 



A. 



(93) 



T] — Tjre.c 



The two new factors in this expression, D{ri) and A, require clarification. 

The factor D[ri) compensates for gradual worsening of our approximation when the wavenumber 
n approaches n^,. For large n's, the functions under the integrals change sign within the window 
^^dec5 instead of being constant there. This leads to the decrease of the true value of the integral 
in comparison with the approximated one. The evaluation of this worsening suggests that it can 
be described by the damping factor 



D{n) 



1 + 



nAr]de 



We inserted this factor 'by hand' in Eq. ()93p. Additional arguments on this damping are given in 
Ref. 3 (see also H)- 

The factor A is the result of the remaining integration over r] in Eq. (El |3> 




dr]'e~T5-M^ \=}^ I dr] 



Vr 
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Since e rapidly changes from to 1 around recombination, the integrand (e" 



-(3/10)r(,7) 



-T{ri) 



) 



is nonzero only there, and A is expected to be of the order of Arjdec- To give a concrete example, 
we approximate g{r]) by a gaussian function 



div) = , — 

V2^(A7?rfec/2) 

Then, the quantity A can be found exactly, 

+00 
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0.96A7]dec- 



Clearly, factor A in Eq. H93|) demonstrates the fact that the CMB polarization is generated only 
during a short interval of time around recombination. 

Finally, substituting (|93|) into l|86jl and recalling l|14j) . we obtain the polarization angular power 
spectrum: 

2 



•' s=l,2 



s 

d h {n,r}) 




drj 


n — ^rec 



PiiCr 



(94) 



Similarly to the case of temperature anisotropies, the projection factors P^(Cdec) predominantly 
translate n into i according to Eq. ()9U() . The oscillatory features of the power spectrum of the first 
time-derivative of metric perturbations get translated into the oscillatory features of the power 
spectra for E and B components of polarization. This is illustrated in Fig.[7|^c) and Fig. [7|^d). The 
waves with nr?^ec did not enter the Hubble radius by ?? = r/^ec- They have no power in the 
spectrum of dh{n,ri)/dr] at n <C 90, and therefore there is no power in polarization at ^ ^ 90. On 
the other hand, the first gravitational peak at n ~ 90 gets reflected in the first polarization peak 
at£^ 90. 

There is certain difference, however, between the projection functions Ei[C,dec) and Bi{C,ciec)- 
This is shown in Fig. El The B-mode projections are more 'smeared' and their maxima are shifted 
to somewhat lower ^'s. This explains the visible difference between Cf^ and Cf^ in Fig. [7| (see also 
The polarization angular power spectra, plotted in Fig. [3 were found from exact numerical 
calculations. 



C. Temperature-Polarization cross correlation 

A very special role belongs to the TE cross-correlation spectrum. We will show below that 
the TE correlation at lower €s must be negative for gravitational waves and positive for density 
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FIG. 6: Relative contributions of an individual Fourier mode n — 100 to various multipoles i in the power 
spectra i{l + l)Cf for temperature and polarization. The normalization has been chosen arbitrarily but 
same in both the graphs. 




FIG. 7: The left panel shows (a) the power spectrum of temperature anisotropics (.{t + l)Cj"^ (in /iK^) 
generated by (b) the power spectrum of g.w. metric perturbations Ijl^l) . (3 = —2. The right panel shows (c) 
the power spectra of polarization anisotropics + \)Cf^ (solid line) and l{t + \)Cf^ (dashed line), panel 
(d) shows the power spectrum of the first time derivative of the same g.w. field, Eq. H15|l . 
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perturbations. This distinctive signature can turn out to be more valuable for identification of 
relic gravitational waves than the presence of the B polarization in the case of gravitational waves 
and its absence in the case of density perturbations. The expected TE signal from gravitational 
waves is about two orders of magnitude stronger than the BB signal, and it is much easier to 
measure. At lower ^'s, the contributions to TE from gravitational waves and density perturbations 
are comparable in absolute value, so the g.w. contribution is not a small effect. The total TE 
cross-correlation has already been measured at some level P, |3| • 

To find the TE power spectrum we have to use the product of aj^{n,s) and af^{n,s) in Eq. 
(jHHJ. For qualitative analysis we will operate with the approximate expressions (|88j) and (|9.S|) . 
Then, the TE correlation reads 

(s s 
A / ^dh{n,r]) s d h* {n,r]) 

(Te{(dec)EiiCdec)) . (95) 

For a given £, the projection factor {T£{Qec)E£{Cdec) ) peaks at n ^ i and is positive there. There- 
fore, the sign of Cf^ is determined by the sign of the term: 



V=Vde 



^ h* t \dh{n,vi) s dh*{n,r]) 



dh?'{n, if) 
dij 



(96) 

V=Vdec: 



V=Vdec 

The adiabatic decrease of the g.w. amplitude upon entering the Hubble radius is preceded by 

s 

the monotonic decrease of the g.w. mode function hn (v) ^ ^ function of ry. This behaviour 
is illustrated in Fig. ^ It is clear from the graph that for n < 100 the quantity (|96|) is negative, 
because the first derivative of hn [f]) is negative. Therefore, for ^ ^ 90 the correlation Cj^ must be 
negative. For larger £'s the TE correlation goes through zero, changes sign and oscillates reflecting 
the oscillations of the function H96|) in the n-space. An exact numerical graph for TE correlation 
caused by gravitational waves is depicted in Fig. |HJ The graph clearly shows how the sign and 
features of the spectrum ()96() get translated into the sign and features of the TE correlation. 

It is shown in Appendix El that in the case of density perturbations the TE correlation must 
be positive at lower €s. This is because the relevant metric perturbations associated with density 
perturbations are growing in time and therefore the first time-derivative of metric perturbations is 

aitive. Because of other contributions, the TE correlation is expected to change sign at £ ~ 70 
, The region of intermediate multipoles 15 < £ < 90 should be of a particular interest. On 
one hand, the multipoles £ > 15 are not affected by the reionization era and its uncertain details. 
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FIG. 8: The bottom panel shows the spectrum H96() of gravitational waves, whereas the top panel shows 
the angular power spectrum £{£ + \)Cj^ caused by these waves. The negative values of these functions are 
depicted by broken lines. 

except the overall suppression by e~^'^''"=»°". On the other hand, at ^ < 90 the g.w. contribution to 
TE is not much smaller numerically than the contribution from density perturbations. The lower 
multipoles ^ ^ 15 are affected by reionization, and we shall study reionization in the next Section. 



VII. EFFECTS OF REIONIZATION ERA 



The reionization of the intergalactic medium by first sources of light has occurred relatively 
late, at z ~ 30 — 7 (see, for example, 0,13]). In contrast to the recombination part of the visibility 
function g{ri), which is narrow and high, the reionization part of g{rj) is broad and much lower 
(see Fig. I15j) . For a crude qualitative analysis, one can still apply analytical formulas derived for 
recombination. One has to replace rjdec with ijreion and Ar/^ec with Aijreion- The waves that start 
entering the Hubble radius at rjreion, i-e. waves with wavenumbers nijreion ~ provide the most 
of power to the spectrum of first time derivative of metric perturbations. Therefore, these waves, 
with n ~ 12, produce a 'bump' in polarization spectra at the projected £ ~ 6. However, these 
wavenumbers are comparable with the wavenumber n**, n^:^:Ar]reion ~ 27r, that characterizes the 
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FIG. 9: The source function ^n{v) rcionization era for different wavenumbers. 

width of the visibihty function at reionization. The assumption n <^ n** is not well satisfied, the 
analytical approximation becomes crude, and one has to rely mostly on numerics for more accurate 
answers. 

The numerical solution to the integral equation (jHJ) in the reionization era is shown in Fig. IHl 
One can see that the damping effect is expected to commence from n > 20, as the polarization 
source function ^niv) begins to show an oscillatory behaviour. The projection relationships are 
also far away from the almost one-to-one correspondence (. k, n that was typical for recombination 
era. In Fig. we show the contributions of a given n to various €s in the polarization power 
spectra + 1)C^. One can see that a considerable portion of power from a given n is distributed 
over many lower-order .^'s. 

The total effect of the reionization era is shown in Fig. 1111 This numerical result was based on 
our simplified model of homogeneous reionization, as described in Appendix^ The reionization 
'bumps' at lower multipoles for E and B components of polarization are similar in shape and 
numerical value. 
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FIG. 10: The contributions to the angular power spectra from an individual mode n. The dashed line is for 
e{£ + l)Cf^ and the solid line is for i{i + l)Cf ^. 
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FIG. 11: The rcionization 'bump' in the polarization power spectra. The dashed lines show the spectra 
when the rcionization era is ignored. The solid lines include both eras. The darker lines are for _B-mode, 
and the lighter lines for i?-mode. 
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VIII. COMPARISON WITH AVAILABLE OBSERVATIONS; SIGNATURES OF RELIC 
GRAVITATIONAL WAVES 



The theory that we are using here is appHcable to any primordial spectral index n. The initial 
conditions for gravitational waves, Eqs. H29|) . H30|) . as well as analogous initial conditions for density 
perturbations, hold /3, and hence n, as a free constant parameter. The spectral index n can be 
larger, equal, or smaller than 1. However, from the theoretical point of view, the 'red' primordial 
spectra n < 1 (/? < —2) seem to be unacceptable, or at least questionable. If /? < —2, the mean 
square fluctuations of gravitational field, Eq. H12|). are power-law divergent in the limit of very 
small wavenumbers n. One could argue that the extrapolation of the primordial spectrum to the 
'infrared' region of very small n's is uncertain, and for some reason the shape of the spectrum 
bends in the 'infrared' region making the integral 1)12(1 convergent at the lower limit. We prefer 
not to hide behind this possibility. If the shape of the primordial spectrum is allowed to be 
varied, then practically anything in the CMB data can be explained by the properly adjusted 
primordial spectrum. Therefore, our theoretical preference (unless the data will enforce us to 
change this preference) is a constant primordial spectral index n > 1 (/3 > —2). Obviously, such 
primordial spectra entail no difficulty in the 'ultraviolet' region of very large n's, because such 
short-wavelength fluctuations (today's wavelength ~ 3 cm.) did not satisfy the requirements of 
superadiabatic amplification and simply have not been generated. 

We now return to the TE correlation. In Fig.Elwe show the contributions to the TE correlation 
function from relic gravitational waves (gw) and primordial density perturbations (dp). (In order 
to enhance lowest ^'s we use the combination {i + l)C£ rather than i{£ + l)Ci.) For this illustration, 
we choose a flat primordial spectrum n = 1 (/? = — 2) and assume equal contributions from (gw) 
and (dp) to the temperature quadrupole: R = 1, where 

R^CfJ,i9w)/CjJ,{dp). 

We include the effects of reionization according to the model with Treion = 0.09. The (gw) contri- 
bution is numerically calculated from the solution to the integral equation ()61() . whereas the (dp) 
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3l|. One can see from the graph that 



contribution is plotted according to the CMBfast code 
the negative TE correlation function at lower ^'s is only possible if there is a significant amount of 
primordial gravitational waves. One can also see from the graph that a mis-interpretation of the to- 
tal TE eff^ect as being caused by density perturbations alone, could lead to a serious mis-estimation 
of 



43 



451 shows clusters of 



The TE correlation at lower £'s measured by the WMAP mission 
data points, including the negative ones, that lie systematically below the theoretical curve based 
on density perturbations alone. It is true that the data points in the interval 10 < ^ < 70 are 
concentrated near a zero level, the error bars are still large, and the measured TE correlation can 
be appreciably different from the theoretical statistically averaged TE correlation. However, the 
recent paper page 26, explicitly emphasizes the detection of the TE anticorrelation by WMAP: 
"The detection of the TE anticorrelation near ^ ~ 30 is a fundamental measurement of the physics 
of the formation of cosmological perturbations...". As we have already stated several times, the 
TE anticorrelation at lower ts, such as ^ ~ 30, can only take place (within the framework of all 
other common assumptions) when a significant amount of relic gravitational waves is present. 

Our theoretical position, as explained in this and previous papers, is such that we are asking 
not 'if relic gravitational waves exist, but 'where' they are hiding in the presently available data. 
We shall now discuss some models that fit the CMB data and contain significant amounts of relic 
gravitational waves. More accurate observations with WMAP, and especially Planck, should firmly 
settle on the issue of the sign and value of the TE correlation at lower multipoles. Hopefully, these 
observations will establish the presence of relic gravitational waves beyond reasonable doubts. 

To sharpen the discussion of allowed parameters, WG take the model with, n — 1.2, '^reion — 

0.09 

and R = l. We take the values ^mh? = 0.03, l^^/i^ = 0.12, h = 0.75, and we normalize the g.w. 
contribution to ^^tii = 440 uK^ at ^ = 2. In Fig. [T3k we show our calculation of the TT 

Q n 

correlation function in comparison with WMAP data and the best fit ACDM model [1^. (If it 
comes to the necessity of explaining 'dark energy', natural modifications of general relativity will 
be superior to unnatural modifications of the matter sector [s^.) One can see from the graph that 
even this model (which lies, arguably, on a somewhat extreme end) is consistent with the TT data 
at all ^'s and significantly alleviates the much discussed tension between theory and experiment at 
^ = 2. In Fig. 113b we show the TE correlation for exactly the same model. One can see that the 
inclusion of relic gravitational waves makes more plausible the negative data points at lower €s. 
Since the relative contribution of gravitational waves becomes small at ^ > 90, the higher £ portion 
of the graph is governed by density perturbations alone. Obviously, models with a little smaller n 
or R reach the same goals. 

Finally, in Fig. we combine together all correlation functions induced by relic gravitational 
waves. The graphs are based on the discussed model with n = 1.2, i? = 1. The future detection 
of the BB correlation will probably be the cleanest proof of presence of relic gravitational waves. 
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FIG. 12: The dotted line shows the contribution of density perturbations alone, and the dashed line shows 
the contribution of gravitational waves alone. The solid line is the sum of these contributions. It is seen 
from the graph that the inclusion of g.w. makes the total curve to be always below the d.p. curve. 
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FIG. 13: The dashed line shows the best fit ACDM model without gravitational waves. The solid line shows 
a model with spectral index n = 1.2 and gravitational waves R = 1. 



0.1 /iK^ in the 



For the discussed model, the predicted level of the BB signal is ^^'^p-Cf^ 

region near the first peak at ^ ~ 90. This level of i3-mode polarization shouiO De aetectaDie Dy 
the experiments currently in preparation, such as CLOVER 51], BICEP 52] and others (see also 
[^). Obviously, in this paper, we ignore many complications, including astrophysical foregrounds. 
More information can be found in Ref. 
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FIG. 14: The summary of CMB temperature and polarization anisotropics due to relic gravitational waves 
with n = 1.2 and R = 1. 

IX. CONCLUSIONS 

In this paper, we summarized the properties of relic gravitational waves that are directly rele- 
vant to CMB calculations. We explained the reasons why we are working with g.w. backgrounds 
possessing specific statistical properties, amplitudes, spectral indices, etc. Then, we worked out, 
essentially from first principles, a theory of CMB anisotropies induced by relic gravitational waves. 
Some parts of this theory are rederivations and confirmations of previous studies, some parts are 
new. The important advantage of our approach, as we see it, is a transparent physical picture. We 
believe we have demonstrated in the paper that every detail of the derived correlation functions is 
under full analytical control. Clear understanding of the participating physical processes has led 
us to the conclusion that the TE correlation in CMB can be a valuable probe of relic gravitational 
waves. We compared our theoretical findings with the WMAP data. We believe that the TE 
anticorrelation detected by WMAP at ^ ~ 30 is certain evidence for relic gravitational waves in 
the already available data. We propose more accurate observations of the TE correlation at lower 
ts and believe that these observations have the potential of providing a firm positive result. 
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APPENDIX A: POLARIZATION STATES AND RANDOMNESS OF GRAVITATIONAL 
WAVES 



As stated above, the quantities | hn (vn describe the magnitude of the mean-square fluctuations 
of the g.w. field in the corresponding polarization states s. We consider a particular mode n of the 
field: 

/iij(n,r/,x) = Pij (n) U„ (??)e*"-'' +hn (f?)e~*"''' c„ j . (Al) 

s=l,2 ^ ^ 

The lowest-order independent g.w. correlation functions amount to 

{hij{n,rj,^)p'^{n) /i^^n', 7?, x') P '^'(n')) = 8 |L {r])f 5,s'S^'\n - n'). (A2) 

The action of the random g.w. field (jAlj) on free particles leads to their relative oscillatory 
motion. We refer this motion to a local inertial frame. Let a ring of free particles to lie in the 
(1, m)-plane; the ring encircles the axis n. Then, the mean-square amplitude of oscillations in the 
'+' polarisation state is determined by | hn P, whereas the mean-square amplitude of oscillations in 

2 

the 'x' polarisation state is determined by | /in p. In general, the random gravitational- wave field 

12 12 

can be such that the oscillation amplitudes | /in P and | /in P are different. But if | /in P / I /in 
then the (averaged) observed picture of oscillations is not symmetric with respect to rotations 
around n. 

Formally, the correlation functions of the field with | /in | 7^ | /in | do not have symmetry with 
respect to the change of polarisation basis. Indeed, the transition to the primed basis according to 
Eq. brings the gravitational wave mode (jAll) to the form 

/i,,(n,r/,x) = Y P'^J (n) (bn (^?)e*"-" + L*(r/)e-*"-^) , 

s=l,2 ^ ^ 

where 

1 1 1 2 2 

bn iv) = cos2ii hn {ri) Cn + sin 2V' hn iv) Cn, 

2 1 1 2 2 

bn iv) = - sin hn (i?) Cn + COS 2ip hn iv) Cn • 
Taking into account the relationships (jlOjl . we can now derive the correlation functions for the 
new polarisation components: 

{bnb*n') = (^\ hn l"^ cos^ 2tP + \ hn P sin^ 2V'^ 6^^^n-n'), 



{bj*n') = (^1 hn psin2 2V' + | L pcos2 2v) S'-^^n-n'), 

1 2 / 1 2 \ 

{bnb*n') = - I /in P - I /in M sin 2-0 COS 2^^ 5(3)(n-n'). 



(A3) 
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It is seen from (IA3|) that, in general, the ■0-dependence survives, and the assumption of statistical 
independence of polarisation components in one basis is not equivalent to this assumption in another 
basis. However, one recovers the original correlation functions <\A2^ from (|A3|) if the conditions 
((TH|) are fulfilled. 

Similar properties hold true for circular polarizations. A g.w. mode /ijj(n, 7y,x) expanded over 
circular polarisation states is given by 

/iij(n,ry,x) = hij{n,r],x) + /i*(n, 77, x), 



where 



hij{n,7],x)= (n) L (??)e*"''' Cn . 



(A4) 



=L,R 



We assume that the complex random coefficients Cy^ (s = L,R), satisfy the statistical conditions 
(dEI). 

The relevant independent correlation functions are calculated to be 

(/ii,(n,?7,x)p*^*(n) /i^,(n',7?,x')p''(n'))=4 L (v)^ ^s.'^^'Hn - n'), {s = L,R). (A5) 



If the observer views the motion of test particles from the — n direction, i.e. against the direction 

R ^ 
hn (r/) 



of the incoming gravitational wave, the function 



is responsible for the mean-square am- 



L 



hn iv) 



IS 



plitude of the right-handed (clockwise) rotations of individual particles. The function 

responsible for the left-handed (anti-clockwise) rotations. (For more details about the motion of 

free particles in the field of gravitational waves, see ^^■) 

Expansion HA4() preserves its form under transformations if one makes the replacements: 

h 'n =hn e^'^^ , h 'n =hn c"*^'^. Therefore, the correlation functions HA5|) do not change, regardless 

2 



the value of the amplitudes 



hn iv) — L,R). On the other hand, discrete transformations H22|l 



L R R L 

generate the replacements: h n =hn, h n =hn- Therefore, the sense of correlation functions (|A5|1 
changes from L to R and viceversa. The symmetry between left and right is violated, unless the 
the conditions (1231) are fulfilled. 



APPENDIX B: ASTROPHYSICAL PREREQUISITES 

(i) Ionization History 

The ionization history of the Universe enters our equations through the density of free electrons 
Ne{r]), i.e. electrons available for Thompson scattering. Specifically, we operate with the quantity 
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q{7]): 

g(r?) = aTa{r])Neirj)- (Bl) 
The optical depth r between some instant of time rj' and a later instant rj is defined by the integral 

V 

r(r?,r?')= / dr,"q{v"). 



rj' 



The optical depth from some rj to the present time r]ji is denoted r(r/) and is given by 

Vr 

Tiv) = t{vr,v) = / dr]'q{7]'). 



V 

It follows from the above definitions that r(r/, rj') = T{rj') — T(rj). 

A key role in our discussion is played by the quantity g{rj) called the visibility function: 

9{V) = q{ri)e-^^^'^ = -e'^^'^ , / g{ri)drj = 1 
dv Jo 

The state of ionization is determined by microphysical processes during all the evolution of 



the Universe 



5£, 



59|. For our purposes it is sufficient to focus on two eras: early recombination 
and late reionization. The recombination of primordial plasma into atomic hydrogen and helium 
is accompanied by decoupling of CMB radiation from the rest of matter. (For simplicity, we do 
not distinguish here the notions of recombination and decoupling.) This relatively quick process 
has happened at redshifts around z^^c ~ 1100. Much later, at redshifts around Zreion ~ 

10 9 the 

intergalactic medium has become ionized again, presumably, by ionizing radiation of first condensed 
objects. 



The density of free electrons is modelled 

Ne{ri) = (l 



60| by the expression 



2 J nip \ a{ri) 

where Yp ~ 0.23 is the primordial helium mass fraction, Xe{i]) is the fraction of ionized electrons, 
0,b is the baryon content, and nip is the mass of a proton. In the framework of linear perturbation 
theory it is sufficient to regard the electron density as homogeneous, i.e. depending only on ry. For 
Xe(?y) we use the fitting formula jo^: 

Yr,\~^ / Co \ / nir, \ „^., _1 / Z \C2-1 / a' 



^'^ V 2y \WOOJ \2aTlHPcJ ^ Viooo/ V« ^ ^ ' ^ ' 

where ci = 0.43, C2 = 16 + 1.8lnO,B, and z is a redshift. 
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FIG. 15: The graphs (a), (b), (c) show the quantities q{ri), T{ri), g{r}). The graph (d) is a zooming oi g{rj) 
in the region of recombination. The dashed hne shows a model without reionization. The solid line takes 
rcionization into account, with Treion = 0.1. 

As for the reionization, we assume that it was practically instantaneous and happened at a 
redshift Zj-eion ~ 16. The function X(. is determined by Eq. (|B2|) for z > Zreion and = 1 for 

Z ^ Zreiorf 

To plot the graphs for q{ri), t{i]) and g{ri) in Fig. 1151 we use the above-mentioned parameters 
and rjf, = 0.046 P, It is seen from Fig. that the visibility function g{7]) is sharply peaked 
at the era of recombination. The peak can be characterized by the position of its maximum r]dec, 
and the characteristic width Arjdec- A similar, but less pronounced, peak is present also at the era 
of reionization. 

(ii) Frequency dependence of the Stokes parameters 

As is seen from Eq. ()46() . the frequency dependence of both, temperature and polarization, 
anisotropies is governed by the function f{i)) 3i]. We assume that the unperturbed radiation field 
has a black-body spectrum. 

Ml") = 7 — ^— , 

exp I hu I kT I — 1 
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where T = Ta{rj) and the present-day value of T is T{r]ff) ~ 2.73K. 

It follows from these formulae that /(P) is approximately 1 in the Reyleigh- Jeans part of the 
spectrum, and /{u) varies as hu/kT in the Wein part. In practice, we are usually interested in 
the total intensity, i.e. Stokes parameters integrated over all photon frequencies. This integration 
produces the factor 7: 



7 



dv -^nQ{u)f{u) = -4/0, 



(B3) 



where /q = / ^^^n^iy) is the total intensity of the unperturbed radiation field. The factor 7 
often appears in the main text of the paper. 

APPENDIX C: TWO ESSENTIAL VARIABLES FOR TEMPERATURE AND POLAR- 
IZATION 



It is seen from Eq. H48|) that the frequency dependence of n'-"'"^ is determined by the factor 
/(t')no(t'). Therefore, we can single out this factor and write 

n^^^(??,z^,/i,(/>) = ^/(z>)no(i>)n(^)(??,/x,(/)). 

The equation that follows from 1)48(1 and (|49() reads 
d 



dr] Att 



+ 1 2tt 

^"^^ f f'dll'dcP'PilI,cp-f,',cP')h(^\7^,f,',<P'). 



(CI) 



-1 



The quantity n^-*-) (rj, (p) as a function of 4> can be expanded in a series 



+00 



m=—oo 



The explicit structure of the Chandrasekhar matrix P(^, /u', is given by the expression 
(equation(220) on p. 42): 

P(^, /i', 0') = Q [PW + (1 - ^2)^ (1 _ ^'2)ip(l) + p(2) 

where matrices Q, p(°)(/i; /i'), P^^\f^,4>; fi', (/)'), P(^)(^, /i', 0') read 

/l o\ 



(C2) 

Q 



Q 



1 
y 2y 
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p{0) 



/2 



1 




^ Afifi' cos I 



p(i) 



sm I 



— 2/ism (</)' — (f 



COS {(j)' - (/)) 



J 



p(2) = 2. 
4 



^ /iV^ COS 2 ((/.'- 0) (,Qs 2 - 0) -^2^' sin 2 ((/.'- ^ 

-^'2 cos 2 (0' - 0) cos 2 (0' - 0) ^' sin 2 {(j)' - 4>) 
y /i;u'2 sin 2 (0' — (/)) —fi sin 2 — 0) ^/x' cos 2 — (p) J 



The structure of the Chandrasekhar matrix P(^, (/>; /i', 0') is such that it does not mix the 
m(/)-dependence, that is, 

+ 1 2tt 



1 



I\„im4> ^imd 
~ e 



-1 



Moreover, this integral vanishes for all \m\ > 2. This means that Eq. ((HH) is a homogeneous 
differential equation for all with m ^ ±2. Assuming zero initial conditions at some initial rj, 
we obtain fi^^ = for all m / ±2. Hence, we are left with three functions of {r],iJ,): 

nW(r?,/i,<A)=nW(r?,/i)e±2^f (C3) 

We are now able to show that only two out of the three functions h^^\r],fj,) are independent. 
Indeed, using (|C3|) in Eq. (|C1() we arrive at a system of three linear equations for the components 
of n(-'-)(r/, /i): 



^ + qiv) + 



ni(?7,/x) 



\ 



1 



3 , ^ 



( ^2 \ 



J(7?), (C4) 



Y ±2i/i y 

where T{r\) is the remaining integral that was left over from the last term in Eq. HClf) : 



'^W) = o / t^At' ^lij], - h2{ri, fi') ± ifi'nsir], fi') 



(C5) 
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Despite the complicated appearance, only two of the three equations HC4|) are really coupled. 
Indeed, making linear combinations of equations HC4() it is easy to show that the combination 
2ifi{hi — n2) =F (1 + A'^)^3 satisfies a homogeneous differential equation. Assuming zero initial 
conditions, we derive 

2ifi{hi - 712) =F (1 + = 0. 



As two independent and essential variables we choose 



In terms of a and /? equations of radiative transfer reduce to Eq. (|51j) and Eq. I|52|). and the 
definition (|(]5|) for I{r]) takes the form of Eq. (|53|) . 

Although we have considered a gravitational wave perturbation, the existence of only two es- 
sential variables is a general statement and it applies to density and rotational perturbations as 
well. In general, the perturbed radiative transfer equation contains an arbitrary function /(r/, /i, (p) 
in front of u, rather than a specific combination (1 — ]i'^)e^'^^'^dhn/ drj quoted in Eq. HC1|1 . Function 

/(ry, /X, 4>) can be expanded in a series similar to Eq. (|C2|1 : 

+00 

m=—oo 

Since the scattering integral (second term in the r.h.s of Eq. 1)011) ) vanishes for all \m\ > 2, the 
functions hm\'r],fJ') are fully determined by fmiv^t^) > 2) and describe the temperature 

variations only. To discuss polarization, we have to consider three remaining cases m = 0, ±1, ±2. 
The explicit identification of the two essential variables a{r], fi) and (3{rj, n) for the case m = ±2 

L 

has been given above. Specifically for gravitational waves, /2 = (1 — l^ )dhn f-2 = 

R 

(1 - fi^)d hn {r])/dr]. 

The identification of the two essential variables for the case m = proceeds in a similar manner. 
Having calculated the scattering integral for m = one can show that the equation, analogous to 
Eq. (|C4() . will now read 



d_ 
where 



\ 



n2(??,//) 



h{Ti,ij) + <i{v)io\'n) 



1 



3 / X 



3;u2 - 2 
1 




\ 



+1 
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+1 



Jiv) = ^ I dfJ-' - 2)ni(r/,/i') +ri2(r?,/x') 



-1 

Obviously, ns satisfies a homogeneous equation and can be put to zero. Specifically for density 
perturbations (see Appendix function /o(r/,/i) consists of terms representing gravitational field 
perturbations and the Doppler term arising due to the baryon velocity: 

We have checked that in the case m = ±1 the problem also reduces to only two essential 
variables. The combination ni — n2 =F ifJ-ns vanishes at the zero initial data. 

APPENDIX D: TEMPERATURE AND POLARIZATION ANISOTROPIES CAUSED 
BY DENSITY PERTURBATIONS 

To discover relic gravitational waves in the CMB data we have to distinguish their effects 
from the effects of density perturbations. The theory of temperature and polarization anisotropies 
caused by primordial density perturbations is very much similar to the theory of relic gravitational 
waves. We start from the metric Fourier expansion (|10)) with the polarisation tensors Q. Having 
derived and solved integral equations of radiative transfer in the presence of density perturbations, 
we arrive at our final goal of distinguishing the TE cross-correlations. 

(i) Radiative transfer equations 

The equations of radiative transfer in the presence of a single mode n of density perturbations 
are similar to Eq. ()48() and read 



9 / N . 
^ + qiv) + Ui 
or) 



-(1)/ - i 



/(z>)no(i>) 



(n) — -j^ g(r?)e Vi 



Qiv) 
47r 



(Dl) 



where the extra term e^Vi takes care of the movement of scattering electrons with respect to the 
chosen synchronous coordinate system 

For technical reasons, it is convenient to work with the 'scalar' h{ri) and 'longitudinal' hi{ri) 
polarization mode functions, instead of the original hn {fj) 3l- The relationship between them is 



hn iv) 



l{h{rj)-hi{r,)), hniv) 
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where the wavenumber index n on h{r]) and is imphcit. Both polarization components of 

metric perturbations participate in Eq. (|D1|) . In the frame associated with the density wave, i.e. 

■ ■ * __ 

for n/n = (0,0, 1), the structures e^e^ Pij (see Eq. lOJ) and e'^Vi = —ifJ^vi, depend on /x = cos 6, but 
not on the azimuthal angle <j). 

By argumentation similar to that in Appendix O one can show that a solution (not vanishing 
on zero initial data) to Eq. (jPlj) must have the form 

"n(?y,^) 



f{i')no{v) 



1 



V / 



(D2) 



Substituting Eq. ()D2|) into Eq. HDljl we arrive at a system of coupled equations for a and /5 



d 



drj 



a„(?7,/i) 
/?n(??,At) 



1 /d/i 



2 Vd?/ 
1 



1 



g(r/) 1- P2{^i)]l2 



(D3a) 
(D3b) 



where 



+1 



^liv) = ^ I dfi an{r],fj,), 



-1 
+1 



^2{v) = 2 I ^1-^ (1 -^2(/^)) /?«(?/, At) - P2{fJ-)an{ri,fi) 



(D4a) 



(D4b) 



The quantity Ti is the monopole component of the perturbed radiation field, whereas T2 is the 
quadrupole {£ = 2) component, responsible for the generation of polarization. 

To make contact with previous work we note that the variables a and f3 are closely related to 
the variables At and Ap from Ref. [s]]. Assuming a black body unperturbed radiation field, we 

d. We also note that the mode functions h and hi are related to the mode 

nn 

,|6J| hy h = —2r], /i; = — (h + 67]). Keeping in mind the difference in 
notations, one can verify that equations (|D3a|) and HD3b|) are equivalent to equations (11) in Ref. 

Q. 

(ii) Integral equations and their solutions 



have At = —a, Ap - 
functions h and r] used in 
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A formal solution to equations (|D3a|) and (|D3b|) can be written as 
an{'n,lj) = J dr]'e 





n' p-'riV,-n')--i-ntJ.{ri-ri') 



(D5a) 
(D5b) 



Proceeding in a manner similar to that in Sec. IIVI we substitute (|D5a|) and HD5b|) into HD4a|) and 
(|D4b|) . After certain rearrangements we arrive at two coupled integral equations for Zi and X2: 

V 







1 



^i{v')jo{x) + -l2{ri')j2{x 



V 



+ J dr] e 




1 f dh dhi d^jo > , M / '\ - f \ 



(D6a) 



^2{v) = J dr]'q{r]')e 




+ J drj e 




1 f dh , . dh. , , , 

2 W''"'^*7 1 



2x" 



« / ■ / N , ^(^^kW , I IX f ixdj2 

+ q{v )Vb{v )■ 



dx^ 



dx 



(D6b) 



where the argument of spherical Bessel functions je{x) is x = n{r] — r]'). 

Equations (|D6a() and HD6b|) together with the continuity equations for matter perturbations 



and Einstein equations for metric perturbations (see for example 
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6l 



63, 
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, form a closed system 
, the radiative transfer 



of coupled integro-differential equations. In previous treatments 
equations were presented as an infinite series of coupled ordinary differential equations. 

Similarly to what was done in Sec. lIVBl we can analyze Eqs. ()D6a|) . HD6b|) in terms of expansions 
in powers of n, 



k=0 



T2 = J^xfV^ 

k=0 



The kernels of Eqs. ()D6a|) . l|D6b|) are expanded in powers of x. In the long- wavelength approx- 
imation, i.e. for x ^ 1, we can limit ourselves by first non-vanishing terms. The zero order 
approximation brings us to the equations 
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-— + 2q(r] jvbir] )x 



The solution to these equations is given by 



4'\v) = J dv' 





4"/*'"'<'''>-'"''"'-c"' 





(D8a) 



15 



dhi 
drj 



1 + '^^ j d7]"q{r]")e 



(D8b) 



(iii) Multipole coefficients 

We are mostly interested in the present-day values of a and (3 and, hence, we put ?? = r//j in 
HD5a|) and ()D5b|) . Irrespective of approximation in which the functions ^2{fl) ^^'^ known, 

the multipole coefficients for the radiation field can be found in a way similar to that in Sec. 



«L(^) = {-iYl^/^^{2l + l)6^oaJ{n), 
^fmi^) = (-i)'7\/4^(2^ + l)'5™oaf(n). 



af^{n) = 0, 



where aj{n), af{n) are given by 

Idh 1 dhi d^ 



VR 



aj{n) 



afin) 



drj 



2)! 



1 »?« 



(D9a) 
(D9b) 



and C = n(77B — r?). As expected, in the case of density perturbations all = for m 7^ and 
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aL(n) = 



l32l.l67j. The formal reason for this is that all quantities in HD2j) do not depend 
on (f) (A more detailed exposition can be found in [^, Q). 

Let us start from the temperature multipoles. For illustration, we consider an instanteneous 
recombination. We replace e""^ with the step function, g{rj) with the (5-function, and we neglect 



the I2 term. Then, we get 



aj{n) 



2^1 + ^^6^ ) hiO 



Vr 



+ 



V=Vdec 



drj 



fdh dhi d^ 
\dri dr] dC,"^ 



Vde 
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This expression can be further simphfied by taking the remaining integral by parts. After some 
rearrangements we arrive at the final expression 



aj{n) 



+ 



2v? drf' 6re drj 

1 d'^hi 



2n^ drj^ 



+ Xi + 



1 dhi 
2n dr] 



+ Vb 



dC 



MO 



V=11dec 



VR 



drj 



drj drj^ ' "^^ 



(DIO) 



ride 



For the growing mode of density perturbations in the matter-dominated era, the integrand of 
the remaining integral vanishes. If, in addition, the intrinsic temperature perturbation Ii and 
the plasma velocity Vh are zero at r? = rjdec-, we recover from Eq. (|Dlfl|) the four terms of the full 
Sachs- Wolfe formula (43) in Ref. _68] . We would like to note in passing that the origin of the often 
used, including this paper, combination + l)Ci is in fact a historical accident, when a wrong 
motivation leads to a convenient notation. This combination arises in the essentially incorrect 
formula + 1)C^ = const that can be derived after the unjustified neglect of three out of four 
terms in the original Sachs- Wolfe formula (43) in Ref. ^^I- This clarification is important not only 
by itself, but also for correct physical interpretation of reasons for the rise toward the 'first peak' 
at ^ ~ 220 of the actually observed function + 1)C^. For more details, see jod ]. 

We now turn to polarization anisotropies. We use the zero order approximation ()D8b|l in HD9b|) . 

/ 



1 ({l + 2)\ 



20 \{l-2)\ 



1 -riR 
2 



dr] ge^^"^ / dr]'e 



\ 



dhi 
drj' 



j + 2n dr]"q{r]")e 



-r{v',v") 



Vbiri") 



As g{r]) is a narrow function, the integrand is localized near r] = %ec- By the same argumentation 
that was used in Sec. lVI Bl the above expression can be reduced to quantities evaluated at r] = rjdec-, 



af{n) 



{£ + 2)1 
20 V(^-2)! 



MO 
C2 



dhi , r\ 
dr] J 



(Dll) 



V=Vdec 



where the factors A and A are of the order of the width of g{rj) in the recombination era. Explicitly, 
they are given by 



vr 



VR 



A 



dll giv) / dr]' e w"^ 



A 



dr] g{r]) j drj' e-To^('''^') j drj" g(7?")e-^(^''^"). 


The angular power spectra in the case of density perturbations are described by formulas anal- 
ogous to ()86() . By the steps similar to those in Section fV CI it can be shown that the power spectra 
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are given by 

(jxx' ^ ^ / ndnaf{n)af'{n), (D12) 



vr 

where C for density perturbations is C = \/247rZp/. 
(iiii) Temperature-polarization correlations 

We focus on the TE correlation at the relatively small multipoles, H. < 70. The dominant 
contribution to these multipoles comes from density perturbations that did not enter the Hubble 
radius at recombination, n < 70. It is sufficient to consider the early time evolution of these 
perturbations in the matter-dominated era. 

Restricting our analysis to the growing mode, we write for the metric and matter perturbations 
(see, for example, S^): 

ST 2 2 
h{r]) = Bn = const, hi{ri) = Bn{nrjf, z[^^ = — = - -Bn{nr]f , Vb = - -Bn{nr]f . 



Substituting this solution into (|D10|) . (|D1H1 and taking into account only the lowest-order terms 
in n (the long-wavelength approximation) we arrive at 

aj ~ -Bn MO , (D13a) 

V=Vdec 



Finally, from Eq. (|D12() we obtain for the TE correlation: 



(D13b) 



^ 10 V(^-2)!y 7 " V C 



(D14) 



V=Vdec 

All terms in the above expression are strictly positive. Therefore, the TE correlation caused 
by density perturbations, in contrast to the case of gravitational waves, must be positive at lower 
multipoles ^ ^ 50 where our approximations are still valid. As mentioned before, this difference 
between gravitational waves and density perturbations boils down to the difference in the sign of 
first time-derivative of the associated metric perturbations. 
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